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ABSTRACT 

Recent detections of the starburst galaxies M82 and NGC 253 by gamma-ray telescopes suggest that galaxies 
rapidly forming massive stars are more luminous at gamma-ray energies compared to their quiescent relatives. 
Building upon those results, we examine a sample of 69 dwarf, spiral, and luminous and ultraluminous infrared 
galaxies at photon energies 0.1-100 GeV using 3 years of data collected by the Large Area Telescope (LAT) 
on the Fermi Gamma-ray Space Telescope {Fermi). Measured fluxes from significantly detected sources and 
flux upper limits for the remaining galaxies are used to explore the physics of cosmic rays in galaxies. We find 
further evidence for quasi-linear scaling relations between gamma-ray luminosity and both radio continuum 
luminosity and total infrared luminosity which apply both to quiescent galaxies of the Local Group and low- 
redshift starburst galaxies (conservative P-values < 0.05 accounting for statistical and systematic uncertainties). 
The normalizations of these scaling relations correspond to luminosity ratios of log(Lo i_iooGev/i'i.4GHz) = 

1 .7 ± 0. 1 (statistical) ± 0.2(disper.sion) and log(Lo. 1-100 Gev/i'8-1000 pm) = "4.3 ± 0. 1 (statistical) ± 0.2(dispersion) for a galaxy 

with a star formation rate of 1 Mq yr"', assuming a Chabrier initial mass function. Using the relationship 
between infrared luminosity and gamma-ray luminosity, the collective intensity of unresolved star-forming 
galaxies at redshifts < z < 2.5 above 0.1 GeV is estimated to be 0.4-2.4 xlO"*" ph cm"^ s"' sr"' (4-23% 
of the intensity of the isotropic diffuse component measured with the LAT). We anticipate that ^ 10 galaxies 
could be detected by their cosmic -ray induced gamma-ray emission during a 10-year Fermi mission. 
Subject headings: cosmic rays — Galaxies; starburst — Gamma rays; galaxies — Gamma rays; diffuse back- 
ground 
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1. INTRODUCTION 

Global emission of most galaxies across much of the elec- 
tromagnetic spectrum from radio to gamma-ray energies is 
related to the formation and destruction of massive stars. The 
exceptions are certain galaxies hosting active galactic nuclei 
(AGN). Throughout this work, we refer to galaxies in which 
non-thermal emission is mainly diffuse in origin rather than 
powered by a supermassive black hole as 'star-forming galax- 
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O and B stars radiate the majority of their bolometric lu- 
minosity at ultraviolet wavelengths, which is efficiently ab- 
sorbed by interstellar du st and re-emitted in the infrared (IR, 
iScoviUe & Youni [19831) . The IR emission dominates the 
spectral ener gy distribution of actively star-forming galaxies 



(reviewed bv lSanders & Mirabellll996[) and is a r obust indi- 
cator of star-formation rates (e.g. lKennicuttill998al) . Starburst 
galaxies are distinguished by kpc-scale regions of intense star- 
forming activity and are often identified by their enormous 
IR luminosities, many being classified as luminous infrared 

Lq) and ultra-luminous in- 
) using the def- 



galaxies (LIRGs, Lg-iooo /.tm > 10 



galaxies jui^ mus, z.R-innn um > lO'^^g 
1 proposed bv lSanders & Mirabell (11996b 



frared galaxies (UL IRGs, Ls- 
inition ] 

Massive stars responsible for the IR emission end their lives 
as core-collapse supemovae, whose remnants are believed to 
be the main cosmic-ray (CR, including all species) accelera- 
tors on galactic scales. The energy budget of CRs observed in 
the solar neighborhood is dominated by relativistic protons. 
If ^ 10% of the mechanical energy of the outgoing shocks 
can be transferred into CR acceleration, supernova remnants 
(SNRs) would be capable of sustaining the locally measured 
CR flux ( iGinzburg & Svrovatskiilll964 . 

In star-forming galaxies, CR electrons and positrons (here- 
after simply electrons) spiraling in interstellar magnetic fields 
radiate diffuse synchrotron emission in the radio continuum 
(RC) at a lev el closely related to the total IR luminosity of 



the galaxies dvan der Kruit 
iHelouet al.1 [19851: iCondoii 



1971 



[1971 [de Jong et al.l[T98l 
19921) . Thermal bremsstrahlung 



emission originating from H II regions ionized by the mas- 
sive stars contributes to the RC emission to a lesser extent. 
Remarkably, the nearly linear empirical RC-IR correlation 
spans over 5 orders of magnitude in luminosity (Condon et al] 
[1991 . 

CRs also create diffuse gamma-ray emission. The same CR 
electrons responsible for radio synchroton emission can pro- 
duce high energy radiation either through interactions with 
gas (bremsstrahlung) or interstellar radiation fields (inverse 
Compton scattering). Inelastic collisions between CR nu- 
clei and ambient gas lead to the production of gamma-rays 
through tt" decay, and also to the production of secondary CR 
leptons by tt^ decay. It is therefore natural to look for correla- 
tions of gamma rays with the CR-induced emissions at lower 
frequencies (and accordingly, IR emissions). 

Two of the nearest starburst galaxies, M82 and NGC 253, 
have been d etected in high-en ergy gamma rays by both 



space-based (lAbdo et alj 



telescopes ( lAcciari et al 



2010a[) and imaging air-Ch erenkov 
TZOOi lAcero et al.l [2009[) . Al- 
though all galaxies are expected to produce CR-induced 
emission at some level, large numbers of SNRs to- 
gether with dense interstellar gas (average number densi- 
ties ^500 cm"^) and intense radiation fields led several au- 
thors to anticipate the central starbursts of M82 and NGC 
253 as detectable gamma-ray sources (e.g. Paglione et alj 
1996; Blom et al. 1999; Domingo-Santamaria & Torres 20051 
Persic et al .1120081: Ide Cea del Pozo et al.l2009tlRephaeh et all 



2O10ir Because the massive stars (M > 8M0) which ulti- 
mately result in core-collapse supernovae have lifetimes of 
0(10 yrs), and particle acceleration in the vicinity of a super- 
nova happens for a short time 0(10"* yrs), the number of CR 
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accelerators in a given galaxy is thought to be closely related 
to the contemporaneous star-formation rate (SFR). Indeed, the 
observed gamma-ray fluxes from M82 and NGC 253 imply 
enhanced CR energy densitities within galaxies undergoing 
rapid massive star-formation. 

A study of the Local Group galaxies detected at GeV 
energies, including the MiUcy Way (MW), Small (SMC, 



Abdo et al 



Abdo et al 



1201 Obi) and Large Magellanic Clouds (LMC, 
120 lOd) . and M31, identified a simple power law 



relation between star-formation and gamma-ray luminosity 
(|Abdo et al. 2010d). On spatial scales within an individual 
galaxy, resolved images of the LMC at GeV energies show 
that gamma-ray emissivity per hydrogen atom, and hence CR 
intensity, is greatest near the massive star-forming region 30 

Doradus ( lAbdo et al.ll2010d). 

EGRET observations (iCillis et al.ll2005h yielded flux upper 
limits for a collection of star-forming galaxies beyond the Lo- 
cal Group (typical limits of 3-5 x 10"** ph cm"- s~' in the 
> 0.1 GeV energy range), and stacking searches for a collec- 
tive sig nal from the sa me galaxies produced no significant de- 
tection. iLenain & Wa lter (2011) recently reported flux upper 
limits in the 0.2-200 GeV energy range for several galaxies 
located within 5 Mpc including M81, M83, IC 342, Maffei 
1, Maffei 2, and M94 using 29 months of data collected by 
the Large Area Telescope (LAT) on board the Fermi Gamma- 
ray Space Telescope {Fermi). The MAGIC Collaboration has 
reported a flux upper limit in the energy range >160 GeV 
towar ds the nearest ULIRG, Arp 220 (15 hours. [Albert et al.l 
120071) . and H.E.S.S. observations of NGC 1068 pro duced a 
flux u pper limit above 210 GeV (4.3 hours. lAharonian et al.l 
I2OOI . 

In this paper, we use three years of Fermi LAT data to per- 
form a systematic search for high energy gamma-ray emission 
from 64 star-forming galaxies beyond the Local Group se- 
lected on the basis of their present star-forming activity. The 
next section describes our galaxy sample. Section |3]outlines 
our analysis of LAT data and we present results in Section 
m including updated spectral energy distributions for signifi- 
cantly detected galaxies and flux upper limits for the remain- 
ing candidates. Five Local Group galaxies previously stud- 
ied in LAT data, the SMC, LMC, Milky Way, M31, and M33, 
are additionally included in a population study of low-redshift 
star-forming galaxies. Our primary objective is to explore the 
global properties of galaxies related to their CR-induced emis- 
sions. We find that a simple power law relationship between 
gamma-r ay luminosity and SFR reported for Local Group 
galaxies (I Abdo et al.ll2010dl) also describes the larger set of 
star-forming galaxies examined here. The implications of this 
scaling relation both for the physics of CRs and the contri- 
bution of non-AGN-dominated star-forming galaxies to the 
isotropic diffuse gamma-ray background are discussed in Sec- 
tion |5] and we predict which galaxies might be detected over 
the course of a 10-year Fermi mission. A standard ACDM 
cosmology with = 0.3, J^a = 0.7, and Ho = 75 km s"' 
Mpc"' is used throughout. 

2. GALAXY SAMPLE 

Massive star-formation is fueled by the availablity of dense 
molecular gas in the interstellar medium (ISM). In order to se- 
lect a sample of galaxies with unambiguous ongoing star for- 
mation, we base our sam ple of galaxies on the HCN survey of 
lGao&Soloinoril ( l2004al) . HCN 7 =1-0 line emission is stim- 
ulated in the presence of dense molecular gas (hh, > 3 x 10"* 



cm ^) typically associated with giant molecular clouds wher e 
the majority of star formation occurs (ISolomon et al.|[T979l) . 
The total quantity of molecular gas in the dense phase as 
measured by HCN line emission exhibits a low-scatter lin- 
ear correlation with both tot al IR lum inosity ( Gao & SolomonI 
l2004bi) and RC luminosity (iLiu & G ao 2010) and is consid- 
ered a reliable indicator of SFR. Observations of Galactic 
molecular clouds suggest that these relations hold at the sca le 
of individual molecular cloud cores as well (IWu et al.ll20()5l) . 

The HCN survey is the most complete study to date in terms 
of total dense molecular gas content, including nearly all of 
the nearby IR-bright galaxies with strong CO emission in the 
northern sky (S > -35°)|3 Additional galaxies with glob- 
ally measured HCN emission data available in the literature 
are also included in the sample. Objects at Galactic latitudes 
j/?! < 10° are excluded due to the strong diffuse emission from 
the Galactic plane. Our final candidate list of 64 galaxies 
beyond the Local Group includes more than a dozen large 
nearby spiral galaxies, 22 LIRGs, and 9 ULIRGs. Global 
properies of the galaxies including RC luminosity, total IR 
luminosity, and HCN luminosity are provided in Table [T] 

Many of the IR-bright starburst galaxies in our candi- 
date list have been previously suggested as interesting tar- 
gets for gamma-ray telescopes on the basis that a sig- 
nificant fraction of CR protons may intera ct in dense 
molecular clouds before escapin g the ISM dvsig [19891 
Akvuz et a lj 119911: iPaglione et alj Il9"96t iTorres et all l20i 
Thompson et al. 1 120071) . furthermore, galaxies well-known 
to host non-thermal leptonic populations are naturally in- 
cluded in this sample since diffuse radio emission associ- 
ated with synchrotron radiation corre lates with IR luminosity 
(iCondon et al.|[T99llYun et al.llIOOl . 

Some of the galaxies in our sample host radio-quiet AGN 
with low-level jet activity. The galaxies associated with 
sources in the Swift BAT 58-mont h survey catalog (15-195 
keV) with AGN-type designations ( iBaumgartner et alll2010l) 
are listed in the rightmost column of Table [T] Hard X-ray 
emissio n is a strong and rela tively unbiased identifier of AGN 
activity ( iBurlon et al.ll201 fl) . while soft X-ray emission from 
AGN is more often obscured by intervening circumnuclear 
and ISM material. X-ray emission from AGN does not re- 
quire presence of a relativistic jet. 

Five Local Group galaxies previously examined in LAT 
data are includ ed in the multiwavelength comparisons appear- 
ing in Section 14.31 Multiwavelength data for those galaxies 
are summarized in Table |2] 

The largest redshift of any galaxy in our sample is z ~ 0.06. 

3. LAT OBSERVATIONS & DATA ANALYSIS 

The Fermi-LAI is a pair-production telescope with large ef- 
fective area (~^8000 cm^ on axis for E >l GeV) and field of 
view (^2.4 sr at 1 GeV), sensitive to gamma rays in the en- 
ergy range from 20 MeV to > 300 GeV. Full details of the 
instrument and descriptions of the on-boa r d and ground data 
processing are provided in lAtwood et al.l (12009b . and infor- 
mation regarding on-orbit calibration procedures is given by 
lAbdo et al. (2009b). The LAT normally operates in a scan- 
ning 'sky-survey' mode which provides coverage of the full 
sky every two orbits (~3 hours). For operational reasons, 
the standard rocking angle (defined as the angle between the 

The HCN survey includes galaxies with 60 /im/100 fim emission larger 
than 50 Jy/100 Jy and CO line brightness temperatures larger than lOOmK for 
spirals or 20mK for LIRGsAJLIRGs. 
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zenith and center of the LAT field of view) for survey mode 
was increased from 35° to 50° on 3 Sep 2009. 

This work uses data collected in sky-survey mode from 
2008 Aug 4 to 2011 Aug 4 for the analysis of 64 celestial 
15° X 15° regions of interest (Rol) centered on the positions 
of the galaxies in our sample. We accept only low -background 
'source' class photon candidate events (lAtwood et al.l 120091) 
corresponding to the P7V6 instrument response functions 
with reconstructed energies 0.1-100 GeV. In order to reduce 
the effects of gamma rays produced by CR interactions in the 
upper atmosphere ( Abdo et al. 2009a), we discard photons ar- 
riving from zenith angles > 100°, exclude time periods when 
part of the Rol was beyond the zenith angle limit, and also 
exclude time periods when the spacecraft rocking angle ex- 
ceeded 52°. 

The data were analyzed using the LAT Science Tools 
software package (version 09-25-02)0 The model for 
each celestial Rol contains templates for the diffuse Galac- 
tic foreground emission (gal_2 Yearp7v6_v0 . fits), a 
spectrum for the isotropic diffuse emission (composed of 
both photons and residual charged particle background, 
iso_p7v6source.txt), and all individual sou rces re- 
porte d in the LAT 2-year source catalog (2FGL, Nol an et aP 
120121) within 12° of the target galaxies. For the individual 
LAT sources, we assumed spectral models and parameters 
reported in the 2FGL catalog. The candidate sources cor- 
responding to the galaxies in our sample were modeled as 
point sources with power-law spectra, dN/dE cx , at the 
optically-determined positions of the galaxies. Due to their 
typical angular sizes of 0(0. 1 °) or smaller, and fluxes close to 
or below the LAT detection threshold, the galaxies considered 
in this work are not expected to be resolved as spatially ex- 
tend ed beyond the ener gy-dependent LAT point-spread func- 
tion jLande et al .1120 121) . 

The normalization and photon index of each gamma-ray 
source candidate were fitted using a maximum likelihood 
procedure suitable for the analysis of binned photon data 
(gtlike). The photons within each Rol were binned spa- 
tially into 0. 1 °-sized pixels and into 30 energy bins uniformly 
spaced in log-energy (10 bins for each power of 10 in energy). 
During the maximum likelihood fitting, the normalizations of 
the diffuse components were left free. We also fit the nor- 
malizations of all neighboring LAT sources within 4° of the 
target galaxy positions, or located within the 15° x 15° Rol 
and detected with exceptionally high significance in the 2FGL 
catalog 

4. RESULTS 

Results from the analysis of 64 galaxies beyond the Local 
Group in LAT data are presented first, followed by a discus- 
sion of multiwavelength relationships using the full sample of 
69 star-forming galaxies. 

4. 1 . LAT Data Analysis Results 

Table [3] summarizes results from the maximum likelihood 
analyses. Source detection significance is determined using 
the Test-Statistic (TS) value, TS = -2(ln(Lo)-ln(Li)), which 

Infomiation regarding the LAT Science Tools package, dif- 
fuse models, instrument response functions, and public data 
access is available from the Fermi Science Support Center 
lhttp://fe rmi ■ gsf c ■ nasa ■ gov/ssc/t - 

'''' More than 1000 attributed photons in the 2FGL dataset or TS > 500; 
see Section|4T|for the definition of TS value 



compares the likelihoods of models including the galaxy un- 
der consideration (Li) and the nu ll-hypothesis of no gamma- 
ray emission from the galaxy (Ln. lMattox et al.iri996l) . Signif- 
icant (r5'>25) gamma-ray excesses above background were 
detected in directions coinciding with four galaxies in the 
sample: two prototypical starburst galaxies, M82 and NGC 
253, and two starburst galaxies also containing Seyfert 2 nu- 
clei, NGC 1068, and NGC 4945. Each of the four gamma-ray 
sources is associated with the correspo nding galaxy from our 
sample according to the 2FGL catalog (iNolan et al.ll2012b . 

We additionally obtain relatively large TS values in the di- 
rections of NGC 2146 (- 20) and M83 (- 15). These ex- 
cesses do not pass the conventional threshold of TS > 25 re- 
quired to claim a source detection. iLenain & Waited (IMTl) 
noted an excess in the vicinity of M83, but determined that 
the best-fit position was inconsistent with the galaxy location 
and instead proposed an association with the blazar 2E 3100. 
We will return to these two galaxies in Section |6] 

CR-induced gamma-ray emission from galaxies is expected 
to be steady on the timescales of our observations. Vari- 
ability is tested for the four significantly detected sources 
by partitioning the full observation period into 12 time inter- 
vals of ~ 90 days each and performing a separate maximum- 
likelihood fit for each of the shorter time-periods. The re- 
sulting lightcurves are shown in Figure [T] None of the four 
sources show significant changes in flux over the three years 
of LAT observat ions. Thi s finding is consistent with the re- 
sults obtained bv fLenain et a l. (2010) for NGC 1068 and NGC 
4945. 

Spectral energy distributions of M82, NGC 253, NGC 
1068, and NGC 4945 are shown in Figure|2] Each flux mea- 
surement represents a separate maximum-likelihood fit fol- 
lowing the procedures described above, but for 6 logarithmi- 
cally spaced energy bins. Each energy bin is further divided 
into 5 logarithmically-spaced sub-bins of energy for binned 
likelihood analysis to maintain ^10 bins for each power of 
10 in energy. Using the power-law spectral models, the maxi- 
mum likelihood photon index values for the four galaxies are 
in the range 2.1-2.4. For M82 and NGC 253, the extrap- 
olated power-law fits from the LAT energy range naturally 
connect with spectral measurements obtai ned by imaging air- 
Cherenkov telescopes at higher ene rgies (lAcciari et al.ll2009t 
lAcero et al.ll2009l: lOhm et al1l20Tll) . 

Integral flux upper limits in the range 0.1-100 GeV at the 
95% confidence level are provided in Table |3]for galaxies not 
significantly detected by the LAT. Limits are computed using 
the profile likelihood technique, for which the flux of a target 
source is varied over a range while simultaneously maximiz- 
ing the model likelihood with respect to all other parameters. 
One-sided 95% CL upper limits correspond to a decrease in 
profile likelihood, Lp, of 2Aln(Lp) = 2.71 relative to the max- 
imum likelihood. The flux upper limits presented in Table |3] 
are derived using a power-law spectral model with photon in- 
dex r = 2.2, corresponding to the typical index of the four 
LAT-detected starbursts. The upper limits increase by 10- 
20% (less constraining) on average when assuming a photon 
index of 2.3. 

4.2. Gamma Rays from Starbiirsting Seyfert 2 Galaxies 

NGC 1068 and NGC 4945 deserve special attention as 
galaxies with both circumnuclear starbursts and radio-quiet 
obscured AGN. The gamma-ray spectra of NGC 1068 and 
NGC 4945 are similar to those of M82 and NGC 253, and 
none of the four galaxies show indications of gamma-ray vari- 
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ability. iLenain et aT] ( 1201 Oh suggested a model for NGC 1068 
in which gamma rays are produced mainly through inverse 
Compton scattering of IR photons by electrons in the mis- 
aligned jet at distances ~' 0.1 pc from the central engine. In 
this scenario, NGC 1068 is viewed as the first member of a 
new class of gamma-ray sources. 

Studies of X-ray selected Seyfert galaxies using LAT data 
demonstrate that Seyfert galaxies hosting radio-qui et AGN 
are generally gamma-ray quiet as a population (Te ng et all 
1201 It lAckermann et al.ll2012bl) . Excepting NGC 1068 and 
NGC 4945, no other radio-quiet Seyfert galaxies have been 
detected by the LAT. The upper limits in 0.1-100 GeV lu- 
minosity inferred for several other nearby Seyfert galax- 
ies are below the gam ma-ray luminosity of NGC 1068 
( lAckermann et al.ll2012bh . The tightest constraint is found for 
NGC 4151, located at a distance of 11.2 Mpc, with gamma- 
ray luminosity < 2 x 10'*'' erg s"' (a factor ~ 10 below that of 
NGC 1068). 

The relative contributions of CR interactions versus AGN 
activity to the gamma-ray emission of NGC 1068 and NGC 
4945 is not yet definitively established. In the multiwave- 
length analysis which follows, we will consider both the full 
sample of analyzed galaxies, and a subsample with galaxies 
hosting Swift-BAT-detected AGN removed. 

4.3. Multiwavelength Luminosity Comparisons 

Many authors have proposed scaling relationships between 
galaxy star-formation rates (SFRs) and gamma-ray luminosi- 
ties motivated by the connections betwe en interstellar gas, 
star formation, supemovae, a nd CRs (e.g. iPavlidou & Fieldsl 
2002HTorres et a l. 2004; Tho mpson et al.l2007HSteckerf2007l 
Persic & Rephaeliii2010i : iLacki et al.ll201 lh ~We now exam- 
ine relationships between gamma-ray luminosity and several 
photometric tracers of star formation taking advantage of the 
three-year accumulation of LAT data. 

The radiative output of star-forming galaxies across much 
of the electromagnetic spectrum is related to the abundance of 
short-lived massive stars. Many photometric estimators of the 
recent star formation histories of galaxies have been used. To- 
tal IR luminosity 8-1000 yum is one well-establis hed tracer o f 
the SFR for late-type galaxie s (reviewed by Kenni cuttll998al) . 
The conversion proposed bv lKennicutll ( Il998h) . 



SFR 

Mo yr-i 



:el.7x 10' 



-10 



(1) 



assumes that thermal emission of interstellar dust approxi- 
mates a calorimetric measure of radiation produced by young 
(10-100 Myr) stellar populations. All SFRs quoted in this 
work consider the stellar mass range 0.1-100 Mq. The factor 
e depends on the assumed initial mass function (I MF), with 
e = 1 for the lSalpete^ (11955b IMF originally used bv lKennicutt 
(11998b'). We will use e = 0.79 to convert to the iChabrier 
(|2003) IMF, proposed after further studies of star formation 
in the 0. 1-1 Mq mass range. 

Several other photometri c SFR e stimators have been cali- 
brated using the iKennicutll (Il998bl) total IR luminosity rela- 
tion. These include RC luminosity at 1.4 GHz pro duced by 
synchrotron-emitting CR electrons (lYun et al.ll2001l) . 



SFR(M0yr-') = e(5.9±1.8) x 10-22Li.4ghz,(WHz-'), (2) 



-h 



We use the conversion factor proposed by ICrain et all ( 120101) for SFRs 
estimated from total IR luminosity: SFRchabricr = 0.79 SFRsalpcter- 



and HCN 7=1-0 line luminosity indicating the quan- 
tity of dense molecula r gas available to form new stars 
dGao & Solomonl2004bl) . 

SFR(Moyr-') = el.8 x 10-%cN(Kkms-' pc^). (3) 

Although these three SFR estimators are intrinsically linked, 
each explores a different stage of stellar evolution and is sub- 
ject to different astrophysical and observational systematic 
uncertainties. 

Figures |3] and |4] compare the gamma-ray luminosities of 
galaxies in our sample to their differential luminosities at 1 .4 
GHz, and total IR luminosities (8-1000 ^m), respectively. A 
second abscissa axis has been drawn on each figure to indicate 
the estimated SFR corresponding to either RC or total IR lu- 
minosity using equations|2]and[T] The upper panels of Figures 
[3] and U directly compare luminosities between wavebands, 
whereas the lower panels compare luminosity ratios. Taken 
at face value, the two figures show a clear positive correla- 
tion between gamma-ray luminosity and SFR, as has been re- 
ported previously in LAT data (see in this context lAbdo et all 
I2010dl) . However, sample selection effects, and galaxies not 
yet detected in gamma rays must be taken into account to 
properly determine the significance of the apparent correla- 
tions. 

We test the significances of multiwavelength correlations 
usi ng the mo dified Kendall r rank correlation test proposed 
bv lAkritas & Siebert (1996). This method is an example 
of 'survival analysis' techniques, suitable for the analysis of 
partially-censored datasets (i.e. containing a mixture of detec- 
tions and upper limits). The Kendall r coefficient, r G [-1 , 1], 
indicates the degree of positive or negative correlation be- 
tween two quanitites by comparing the ordering of each pair 
of points in the dataset. For example, t = 1 describes a set 
of uncensored data points which are monotonically increas- 
ing. This non-parametric approach makes no assumption re- 
garding the particular mathematical form for the relationship 
between compared quantities. See Appendix lAl for a detailed 
description of the method. 

The multiwavelength relationships are tested in luminosity- 
space because we are primarily interested in the intrinsic 
galaxy properties, and importantly, because multiwavelength 
comparisons in flux-space can either falsely suggest a non- 
existant correlation or obscure a genuine physical correlation 
if the underl ying relationship between intrinsic luminosities 
is non- linear (iFeige lson &^erg|[T983t iKembhavi et"aL|[T986l 
and see Appendix We performed a series of Monte Carlo 
simulations to investigate potential biases resulting from se- 
lection effects for flux-limited samples of objects compared in 
luminosity-space. The generalized Kendall t correlation test 
is found to be robust against the selection effects expected for 
our application. A brief discussion of results from that study 
is presented in AppendixIC] 

Significances of the multiwavelength relationships are com- 
puted by comparing the t correlation coefficient of the ac- 
tual data to the distribution of r correlation coefficients which 
could be obtained under the null hypothesis of independence 
between wavebands. Null hypothesis datasets are generated 
by scrambling derived gamma-ray luminosities among galax- 
ies. Only 'observable' permutations are retained to account 
for the truncation of measurable gamma-ray luminosities im- 
posed by the LAT flux s ensitivity (following the method of 
lEfron & PetrosiarJ (11999b ). Specifically, we require that the 
resultant gamma-ray flux (or upper Umit) of each galaxy in a 
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scrambled data set could have been measurable by the LAT in 
order for the scrambled data set to be included in the null hy- 
pothesis distribution. In practice, each permutation represents 
an exchange of gamma-ray luminosities between two ran- 
domly selected galaxies so that statistically meaningful null 
hypothesis distributions can be efficiently generated. Gamma- 
ray flux sensitivity thresholds are determined by randomly 
sampling from the distribution of gamma-ray flux upper lim- 
its presented in Table [3] for each galaxy in the pair to be ex- 
changed. We perform 1000 such exchanges on the actual data 
before counting permutations towards the null hypothesis dis- 
tribution. 

Inaccuracies in the distance measurements to the galaxies 
are included in our analysis because these errors are prop- 
agated in the transformation of measured fluxes into lumi- 
nosities. Scatter in distance measurements applied to a col- 
lection of objects tends to broaden the distribution of objects 
in luminosity-space, and can therefore induce positive corre- 
lation. For the galaxies in our sample at distances 1 Mpc< 
D <100 Mpc, the typical dista nce measurement prec ision for 
individual galaxies is 10-20% dFreedman et al.l20^ll and ref- 
erences therein). For our particular sample, we can compare 
the dist ances provided in the HCN survey of Gao & Solomo3 
(I2004al) to t hose reported in the I RAS Revised Bright Galax- 
ies Sample jSanders et al.ll20()3h . and find that for galaxies 
beyond the Local Group, the average discrepancy between re- 
ported distances is 10%. We assume that measurement errors 
are normally distributed and centered on the true distances to 
the galaxies. 

Figure |5] shows an example of the correlation significance 
test applied to the relationship between gamma-ray and IR lu- 
minosity using the full sample of star-forming galaxies. The 
distribution of correlation coefficients for 10*" 'observable' 
permutations of the actual data are plotted for 3 levels of dis- 
tance measurement scatter The null hypothesis distributions 
tend towards positive values due to the truncation of measure- 
able gamma-ray luminosities, effectively requiring a larger 
correlation coefficient for the actual data in order to claim 
a significant relationship between wavebands. As the level 
of scatter in distance measurements is increased, the tails of 
null hypothesis distributions extend to larger positive values. 
The probability that a null hypothesis data set would have a 
correlation coefficient larger than that of the actual data cor- 
responds to the P-value measure of correlation significance. 

Given uncertainties for the measured gamma-ray fluxes of 
galaxies in our sample, the correlation coefficient of the ac- 
tual data is more appropriately viewed as a probability den- 
sity as opposed to a single value (the Kendall t statistic does 
not explicitly include measurement errors). To account for 
this uncertainty, we draw a flux value for each LAT-detected 
galaxy from a normal distribution centered on the best-fit flux 
with standard deviation equal to the reported flux uncertainty. 
A correlation coefficient is then computed using the sampled 
gamma-ray fluxes. This process is repeated many times to 
generate a distribution of correlation coefficients represent- 
ing the actual data, as shown by the gray bands in Figure |5] 
Finally, we integrate over the probability density of correla- 
tion coefficients representing the actual data to obtain signifi- 
cances for the multiwavelength relations. 

The results of the correlation tests are summarized in Ta- 
ble |4] Using the full sample of 69 star-forming galaxies, the 
null hypothesis of independence between gamma-ray lumi- 
nosity and either RC or total IR luminosity can be rejected at 



the >99% confidence level allowing for 20% uncertainty on 
the distance measurements to the galaxies. No evidence for a 
correlation between gamma-ray luminosity and HCN line lu- 
minosity is found in this study, an apparent discrepancy con- 
sidering the previously establis hed corre lations between RC, 
IR, and HCN line luminosities ( iLiu & Ga o 2010). However, 
global HCN line luminosity estimates are not available for the 
SMC and LMC, M3 1, and M33, which effectively reduces the 
luminosity range over which the correlation between wave- 
bands can be tested, and excludes three LAT-detected galaxies 
from the analysis. 

We also considered the more conservative case in which 
galaxies hosting ^wZ/f-BAT detected AGN, including NGC 
1068 and NGC 4945 along with seven others, were removed 
from the sample given the potential contributions of the AGN 
to the broad-band emission of those galaxies. In that case, the 
P-values for the correlations between gamma-ray luminosity 
and either RC or IR luminosity are ^ 0.05 allowing for 20% 
uncertainty on the distance measurements. 

Even in the most conservative case described above, the 
data prefer a correlation between gamma-ray luminosity and 
both RC and IR luminosity. Although strong conclusions re- 
garding the significance of the multiwavelength correlations 
cannot be made at the time, at minimum, scaling relations de- 
rived from the current sample have predictive value and offer 
a testable hypothesis for further observations. 

We fit scaling relationships between wavebands using sim- 
ple power law forms. For example. 



, ,'i-0.1-100Gev\ , / i'S-lOOOum 

log Ti — =alog' 



ergs 



lOiOL. 



+ (4) 



parameterizes the relationship between gamma-ray and total 
IR luminosity. 

Two regression methods are employed: the Expectation- 
Maximization (EM) algorithm, and the Buckley-James algo- 
rithm. The EM algorithm is similar to the least-squares fitting 
method, and assumes that the intrinsic residuals of gamma- 
ray luminosity are normally distributed in log-space about the 
regression line for fixed values of either RC or IR luminos- 
ity. Non-detections are incorporated in the regression analy- 
sis by determining the degree to which the upper limits are 
compatible with the assumed dispersion about the regression 
line. The variance of the intrinsic residuals becomes a third 
parameter in the fit. The Buckley-James algorithm is similar 
in approach, but generalizes to cases for which the intrinsic 
distribution of residuals relative to the regression line is not 
normally distributed. Instead, the intrinsic dispersion is esti- 
mated using the Kaplan-Meier di stribution obtained from the 
scatter within the dataset itself. Ilsobe et aTl ( II986I) provide 
details regarding the theory and implementation several com- 
monly used survival analysis methods including the EM and 
Buckley-James regression algorithms. 

Table |5] reports the best-fits of gamma-ray luminosity ver- 
sus either RC luminosity or total IR luminosi ty, which were 
evaluated using the ASURV Rev 1.2 code dLaVallev et alj 
1992)0 The two regression methods yield consistent re- 
sults. For both the RC and total IR tracers of SFR, a nearly 
linear power law index of q;=1.0-1.2 is preferred. The scal- 
ing indices between gamma-ray luminosity and SFR tracers 

ASURV and other tools for the analysis of censored data 

sets are available from the Penn State Center for A strostati sties 

Ihttp : / /www .astrostatistics.psu. edu/ st at codes /sc_censor . html) . 
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found for the larger sample of galaxies is consistent within 
uncertainties with that obta ined considering th e Local Group 
galaxies alone, q;=1.4±0.3 ( lAbdo et alj|2010dl) . In particular, 
the gamma-ray luminosity upper limits found for LIRGs and 
ULIRGs disfavor a strongly non-linear scaling relation. Ta- 
ble |5] also contains regression parameters fitted after remov- 
ing galaxies hosting S'w/ff-BAT detected AGN from the sam- 
ple, which would be appropriate if the broadband emissions of 
those galaxies were heavily influenced by AGN activity. The 
fitted parameters in both cases are consistent within statistical 
uncertainties. 

Note that the intrinsic dispersion values presented in Table 
|5]should be viewed as upper limits because we have not made 
an attempt to account for measurement uncertainties in the 
gamma-ray fluxes. 

The best-fit power laws obtained using the EM algorithm 
are plotted in Figures|3]and|4] The darker shaded regions rep- 
resent uncertainty in the fitted parameters from the regression 
and the lighter shaded regions indicate the addition of intrinsic 
residuals to one standard deviation. None of the limits from 
non-detected galaxies are in strong conflict with the best-fit 
relations. 

Luminosity ratios of the form shown in the lower panels of 
Figures [3] and |4] have a direct physical meaning in terms of 
the energy radiated in different parts of the electromagnetic 
spectrum. For a galaxy with non-thermal emission powered 
mainly by CR interactions, and having a SFR of 1 yr"' 
(similar to the Milky Way), the corresponding luminosity ra- 
tios between wavebands are 



log ( ) = 1.7 ±0.1 (statistical) ±0.2(dispersion), (5) 



A A GHz 



and 



log ( r" ' ] =-4.3 ±0.1(statistical)±0.2(dispersion), (6) 

V i^8-1000^im / 

found using the full sample of galaxies with the EM regres- 
sion algorithm. 

Further gamma-ray observations are required to conclu- 
sively establish the multiwavelength correlations examined in 
this section. Additional data may also show that scaling rela- 
tions beyond simple power law forms are required once far- 
ther and more active sources are either detected or are con- 
strained by more stringent gamma-ray upper limits. 

5. DISCUSSION 

In recent years, increasingly sensitive observations of exter- 
nal galaxies at both GeV and TeV energies have shown that 
starburst galaxies such as M82 and NGC 253 are character- 
ized by harder gamma- ray spectra relative to quiescent galax- 
ies of the Local Group (lAce ro et al.'20Q9l: lAcciari et al.ll2009l 
lAbdo et al.ll2010aHOhm et al..i,201 L) . and that global gamma- 
ray lumin osities of galaxies li kely scale quasi-linearl y with 
SFRs (e.g. lAbdo et alJlMOaHlLenain & Walte3l20TTh . 

Figure |6] compares the spectra of eight star-forming galax- 
ies detected by gamma-ray telescopes. The more luminous 
galaxies also have comparatively harder gamma-ray spectra. 
Whereas the power-law spectral indices of M82 and NGC 253 
are 2.2-2.3 extending to TeV energ ies, the spectra of th e LMC 
and SMC steepen above - 2 GeV ( lAbdo et al.ll2010bllg|) . 

Scaling relations of the type examined in Section 14.31 and 
the observed gamma-ray spectra of star-forming galaxies, 



have implications both for the physics of CRs in the ISM and 
for the contribution of star-forming galaxies to the isotropic 
diffuse gamma-ray emission. In this section, we concentrate 
on the global non-thermal radiation features of star-forming 
galaxies from a population standpoint. 

5.1. Physics of Cosmic Rays in Star-forming Galaxies 

A discussion of scaling relations for diffuse gamma-ray ra- 
diation from star-forming galaxies invites comparison with 
the extensively studied RC-IR correlation. The RC-IR cor- 
relation extends over multiple galaxy types including irreg- 
ulars, spirals, and ellipticals with on-going star formation 
(Ide Jong et alJ [1981 iWunderlich et all [19871 |D resse 3 [ma 



iWrobel & Heeschen|[l988[ Il991h . covers galaxies with mag- 
netic energ y densities likely ranging at least 4 orde rs of 
magnitude (ICondon et al.l 119911 : [Thompson et al J I2006I). and 
applie s to galaxies out to at least 7' ^ \ (lAppleton et all 
|2003)- High-resolution imaging has also shown the corre- 
lati on to hold on ~100 pc scales within individual galax- 
ies ([Marsh & Helou|[T995l ; lPaladino et al.ll2006t [Murphy et alj 

A simple model proposed to explain the RC-IR correlation 
posits that galaxies are ef fective 'calo rimeters' of both UV 
photons and CR electrons (IVolklll989h . In this picture, if the 
ISM is optically thick to UV photons, reprocessed starlight 
emitted in the IR becomes a proxy for the abundance of 
massive stars, and hence SFR. The luminosity of CR elec- 
trons is assumed to be proportional to the SFR in the basic 
calorimeter model. If either synchrotron losses dominate for 
CR electrons, or if the ratio of synchrotron losses to other 
loss-mechanisms is uniform for many galaxy types, then the 
RC luminosity should also increase together with the SFR. 
The latter possibility requires some degree of fine-tuning (but 
see lLacki et al.|[2010l) . Numerical models of CR propagation 
in the Milky Way su ggest that our Gala xy is an effective CR 
electron calorimeter ([Strong et al .112010]) . 

Gamma-ray observations extend our study of CRs to in- 
clude hadronic populations beyond the solar system, and no- 
tably, to measure the large-scale propagated spectra of CR nu- 
clei in the Milky Way and other galaxies. At energies above 
the threshold for pion production, hadronic gamma rays have 
the same spectral index as the underlying CR proton spectrum 
in the thin-target regime relevant for proton-proton interac- 
tions in the ISM. Diffuse pionic gamma rays presumably offer 
direct access to the majority of CR energy content, although 
leptonic diffuse emission and individual sources within galax- 
ies must also contribute to the total emission. The gamma-ray 
spectra of the SMC and LMC are both consistent with models 
in w hich pionic gamma -rays dominate the high energy emis- 
sion ([Abdo et al.l2010b l[g[). Many authors have argued that the 
majority of starburst emission at energies > 0.1 GeV is pro- 
duced by interactions of CR nuclei in the interstellar medium, 
assuming primary CR proton to e lectron ratios similar to val- 
ues found in the Milky Way (e.g.[Paglione et anil996l : iTorresI 
[200llPersic et al.ll2008l: [Tacki et al.|[20l'ol) ! 

For the Milky Way in particular, global CR-induced 
gamma-ray em ission above 0.1 GeV stems mostly from neu- 
tral pion decay (lBloemenll985l:lBertsch et al.l 19931) . with lep- 
tonic proc esses contributing 30^0% in the energy range 0.1- 
100 GeV ([Strong et al.ll201(J) . The diffuse neutral-pion-decay 
component for the Milky Way has spectral index ^ 2 .75 above 
^ 1 GeV, matching the locally measured CR proton spec - 
trum (lSimpsonlll983l: iSanuki et al.|[2000l: lAdriani et al.|[20rTI) . 
CR escape from the Milky Way plays a major role in the 
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formation of the measured CR spectrum, but other galaxies 
with higher concentrations of ambient gas may work differ- 
ently, behaving as calori meters of CR nuclei, as suggested by 
[Thompson et alJ (l2007h . Approximately 10% (possibly up to 
20%) of the global gamma-ray emission of the Milky Way is 
thought to originate from individual sources (Strong 2007i) . 

We now review models for the interactions of CR nu- 
clei to examine physical quantities involved in the scaling 
of hadronic gamm a-ray emission, fo l lowin g the formalism 
of recent works by iPersic & Rephaelil (120101) and lLacki et al.l 
dlOll). Neutral pion decay emission is related both to the 
energy density of CR nuclei. Up, and to the ambient gas den- 
sity, «, which serves as target material for inelastic collisions. 
The total hadronic gamma-ray luminosity (erg s~') can be ex- 
pressed as 



Ej-^n{F)Up(r)dE.ydV, 



(7) 



integrating over gamma-ray energies, E^, and over the whole 
galactic volume, V, which comprises the disk and surrounding 
halo occupied by CRs. The factor q{> E^) denotes gamma- 
ray emissivity ( s~' H-atom~' eV~' cm"*) normalized to the CR 
energy density dPrurv et al.lll994l) . Next, we define an aver- 
age effective density of ambient gas encountered by CR nu- 
clei: 



(neff) ■ 



J^n{F)Up(r)dV 

JyUp{F)dV 



(8) 



Note that («eff) can be substantially lower than the average 
gas density found in the galaxy disk. CR nuclei of the Milky 
Way are thought to spend most of their time in the low-gas- 
density halo, based upon combined measurements of unstable 
isotopes and spallation products (e.g. the B/C ratio) in the 
locall y-observed CR composition (reviewed by IStrong et al.l 
l2007h . Using the above definition for the average effective 
gas density, equation]?] separates to 



UpdV X («eff) X 



dq 



dE^. 



(9) 



The leading factor of equation ]9]represents the total energy 
of CR nuclei, and can be re-written as the product of the total 
luminosity of CR nuclei, Lp, and their characteristic residence 
time in the galactic volume, r^s: 



UpdV: 



(10) 



The trailing factor of equation ]9] depends upon the parent 
spectrum of CR nuclei, for which we assume a power law 
spectral distribution, dNp/dKp oc (Kp/l GeV)"'"p, where Kp 
is the kinetic energy. Gamma-ray emissivities (s~' H-atom~' 
cm^) in the energy range 0.1-100 GeV, 



Qo. 



1-100 GeV ■ 



100 GeV 



0.1 GeV 



E,^dE„ 

'dE^ ^' 



(11) 



are computed using the proton-proton interaction cross sec- 
tion parametrizations of iKamae et al.l (]2006t) . Note that the 
emissivities are normalized to the average energy density of 
CR nuclei {\ < Kp < IQ^ GeV) within the entire galactic 
volume. A nuclear enhancement factor of 1.85 is included 



to account f or heavy nuclei in both CRs and target matter 
( iMoril 120091) . The emissivities for Fp = 2 - 3 range from 
Qo.i-iooGeV = (5.8-7.9) X 10-1^ s-' H-atom-' cm\ 
Substituting equations ]T0] and [TTl into equation ]9] yields 



^0.1-100 GeV,7r» - J^-p '^res ("eff) Qo. 



1-100 GeV- 



(12) 



In the paradigm that SNRs are the pri mary sources of galac- 
tic CRs (iGinzburg & Syrovatskiill 19641) . CR luminosity is the 
product of the supernova rate, Fsn, the kinetic energy released 
per supernova, E^n, and the fraction of kinetic energy going 
into CR nuclei with kinetic energies > 1 GeV (i.e. accelera- 
tion efficiency), r/: 



Ep = Fsn fi'sN rj- 



(13) 



Core-collapse supernova rates can be estimated from SFRs 
given an initial mass function and a minimum stellar mass 
required to produc e a core-collapse supernova, which we take 
to be 8 Mq. For a IChabrie"r] (120031) initial mass function, the 
transformation between SFR in the mass range 0.1-100 Mq 
and core-collapse supernova rate is Fsn '5 = SFR, with = 
83M0. The average kinetic energy released per core-collapse 
supernova is ^sn 10^' erg (IWoosley & Weaveill 19951) . The 
gamma-ray luminosity becomes 



i'0. 1-100 GeV,7r» = * ' SFR EsN V 7"res («eff) Qo.l-lOOGeV- (14) 

Greater physical intuition can be realized by inserting values 
for the physical quantities similar to those of the Milky Way, 



LoA-m Gev.TrO = 9.4 X lO^'^crg s ' 



SFR 



Meyr-i 

V \ f Ties \ / (Weff) 



lO^iergy/ VO.I/ VlOV/ \cm 



(15) 



assuming a power law index for CR nuclei of Fp = 2.75 cor- 
responding to Qo.i-iooGeV =7.8 X 10"'^ s"' H-atom"' cm-'. 

A consistency check for this simple approach can be per- 
formed for the specific case of the Milky Way, for which 
the column density of material tr aversed by CR nuc lei, or 
'grammage,' is well known (e.g. iDogiel et al.l 120021) . The 
mean escape length for CR nuclei in the MiUcy Way is .f = 
Tres nip {lien) c « 12 g cm"^ d Jones et al.ll2001l) (cf. the inter- 
action length is ~55 g cm"^). We may equivalently express 
equationJTslin terms of the escape length as 



A).l-100 GeV,7rO = 6.0 X lO^'^Crg S 

EsN 



lO^'erg 



0.1 



SFR 

Moyr- 



10 g cm 



(16) 



If we take SFR =1.6 Mp,yT~^ correspond ing to a SN rate of 
1.9 (±1.1) per century (iDiehl et al.ll2006i) . then our estimate 
for the pionic gamma-ray luminosity of the MiUcy Way be- 



■^0.1-100 



GeV,7r« = 1-1 X lO^Vg S 



^SN 



lO^iergy VO.l 



(17) 



For comparison, the luminosity estimated from a full numer- 
ical treatment (GALPROP code) incorporating an ensemble 
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of multiwavelength and CR measurements is io.i-iooGeV.7r» = 
(5 - 7) X lOf erg s"', with £sn?7 = (0.3- 1) x 10^° erg 
jStrong et alJlMO) . 

It is commonly assumed that the kinetic energy released 
per supernova, and acceleration efficiency are universal from 
galaxy to galaxy. In that case, the model expressed in equa- 
tions [15] and [16] implies that the hadronic gamma-ray lumi- 
nosity of the interstellar medium scales linearly with the SFR, 
and with the product of residence time and effective ambi- 
ent gas density (or equivalently, grammage). These quanti- 
ties are likely correlated in real galaxies. For example, sev- 
eral authors have argued that SFRs and gas densities are re- 
lated (e.g. via the Schmidt-Kennicutt law) so that a non-linear 
scaling of gamma-ray luminosities with SFRs should be ex- 
pected: i-mnr.p v tt" oc SFR' "*-' ^ (iPersic & RephaelillMTol: 
iFields et al.il201oli . 

In general, the residence time of CR nuclei in a galaxy, Tres, 
can be expressed as 



r-' =r-i+r-i. 

res esc pp 



(18) 



A galaxy becomes a 'calorimeter' of CR nuclei when the 
residence time approximately equals the collisional energy 
loss timescale, Tres ~ Tpp, i.e. energy losses are dominated 
by inelastic collisions with interstellar matter rather than by 
escape of energetic particles. The proton-proton collisional 
energy loss timescale depends primarily on the average gas 
density encountered by CR nuclei, and is nearly indepen- 
dent of energy because the inelastic cross-section, cTpp, is 
only weakly energy-dependent for CRs with kinetic energies 
greater than ~ 1 GeV. The collisi onal energy loss timescale is 
(iMannheim & Schlickeisedl 19941) 



Tpp « (0.65 (neff) CTpp c) ' w 5 X lO^yr 



("eff) 



By contrast, the escape energy loss timescale. 



adv 1 



(19) 



(20) 



may be energy-dependent since the diffusion time, Tdif, de- 
pends upon particle rigidity in interstellar magnetic fields. For 
the Milky Way, the main contribution to the gamma-ray lu- 
minosity comes from photons with energies '^l GeV, corre- 
sponding to CR proton energies of ^10 GeV. The diffusive es- 
cape timescale can be estimated as Tdif ^ zl/6D ~ (1 - 2) x 10^ 
yr for a halo height of Zh = 4 kpc, and taking diffusion coef- 
ficients in the range D = (4-8) x 10^*^ c m- s~' for a particle 
rigidity of ~10 GV ( iPtuskin et al.ll2006h . Given a grammage 
of 12 g cm"^, the average gas density corresponding to an es- 
cape time of 2 X 10^ yr is («eff) ~ 0.3 cm"-' (implying that 
Tpp 1 .5 X 10^ yr). In spite of modeling differences, it is gen- 
erally understood that the diffusive escape timescale is much 
shorter than the co llisional loss timescale for the Milky Way 
jStrong et al.ll2007l and references therein). 

Large uncertainties exist for the escape time due to uncer- 
tainty in diffusion parameters, and because CRs might also 
be advected outwards by bulk motion 'supe rwinds' which 
are observed for many starbursts (iLehnert & H eckman 19961. 
The advective timescale, Tadv, like the collisional energy loss 
timescale, is nearly independent of energy. 

The spectrum of pionic gamma rays will be affected by 
a transition between dominant energy loss mechanisms, be- 



coming softer for galaxies in which energy-dependent diffu- 
sive losses are important. In particular, the harder gamma- 
ray spectra found for starburst galaxies may be understood if 
the energy loss rate due to either proton-proton interactions or 
advection is considerably faster than the diffusion timescale. 
A scaling relation for gamma-ray luminosity should include 
the possibility that residence times for CR nuclei in different 
galaxies may vary substantially. 

In the calorimeteric limit, the residence time and effective 
density of target material are inversely proportional, and ac- 
cordingly, the anticipated scaling of pionic gamma-ray lumi- 
nosity becomes a linear scaling of the SFR: 



--0.1-100GeV.7r''krc 



: 5 X lO^'Vg s 



SFR 



TTT • (21) 



M0yr-i 

_EsN_\ /Jl 
lO^iergJ VO.l 

For the calorimeter case, we assume an underlying spec- 
trum of CR nuclei described by a power law with Fp = 2.2 
(Qo. 1-100 GeV =7.8 X 10"'^ s~' H-atom~' cm"*), representative 
of the LAT-detected starbursts. The gamma-ray luminosity 
expected in the calorimetric limit for CR nuclei is plotted in 
Figures [3] and |4] assuming an average CR luminosity per su- 
pernova of ^SN V = 10^° erg [3 

To determine how well the theoretical calorimetic limit 
above describes nearby starburst galaxies, we can transform 
the observed scaling relation between gamma-ray luminosity 
and total IR luminosity to a scaling relation between gamma- 
ray luminosity and SFR using equation[T] This produces 



.^0.1-100 GeV Iscaling relation 



■■N 



SFR 1 

Moyr-i L7^ 



(22) 



where = 1 .^S^^j] x 10^^'' erg s"' and a = 1 . 16 ± 0.07, fit us- 
ing the EM algorithm for the complete sample of 69 galaxies. 
Substituting best-fit values for the normalization and index, 
and assuming a Chabrier IMF (e = 0.79) yields 



.1—100 GeV I scaling relation 



:(1.3±0.3)x lO^Vgs"' 

SPR \ 1.16±0.07 

Moyr- 



(23) 



The ratio between the gamma-ray luminosity estimated via 
the observed scaling relation with total IR luminosity (equa- 
tion [23]) and predicted in the calorimetric limit for CR nuclei 
(equation ITTTi provides an estimate for the calorimetric effi- 
ciency of CR nuclei: 



^ io.l-100Gev|scaling relation 
^ ^'O.l-lOOGeV.Trolrr.s-Tpp 



:(0.3±0.1) 



SFR 

Moyr-i 



0.16±0.07 



VlO^iergy VO.I 



TTt) ■ (24) 



^^ ILacki et alj )201lh calculated the ratio L>i Gev/^s-iooo ^im = 3.1 X 10~ 
expected in the caloiimetric limit for CR nuclei with Fp = 2.0 and fisN*? = 
10" erg. Substituting the con'esponding gamma-ray emissivity, Q>iGeV = 
1.2 X 10"'* s"' H-atom"' cm^. into equation l21l vields L~.i r.,v /Lg_innn ..^ = 
2.5 X 10^, an agreement to within ~20%. 



10 



The Fermi LAT Collaboration 



This calorimetric efficiency estimate should be viewed as an 
upper limit because leptonic processes and individual sources 
must also contribute to the observed global gamma-ray emis- 
sion of galaxies. We may infer from equation|24]that starburst 
galaxies with SFR ^ 10 Moyr"' have maximum calorimet- 
ric efficiencies of 30-50% whereas dwarf galaxies with SFR 
~ 0. IMQyr"' have lower calorimetric efficiencies of 10-20%, 
when assuming an average CR luminosity per supernova of 
EsN V = 10^*^ erg, using a Chabrier IMF, and attributing all 
observed gamma-rays to hadronic processes. Real galaxies 
are expected to exhibit some dispersion about the trend. The 
estimates pr esented here are in good agreement with the cal- 
culations of iLacki et all jlOlll) for M82 and NGC 253. Our 
own Milky Way is the only galaxy for which global gamma- 
ray luminosity can be compared to direct CR measurements, 
and the CR energetics are in fair agreement with our estimates 
above considering that the estimated gamma-ray luminosity is 
below (but consistent with) the best-fit scaling relation. 

Energy-independent cooling mechanisms for CR nuclei ca- 
pable of preserving the injected spectral shape are preferred 
for the LAT-detected starburst galaxies given their relatively 
hard gamma-ray spectra. An enhancement of interaction effi- 
ciency is suggest ed b y the slightly non-linear scaling relations 
found in Section 143] However, the overall normalizations of 
these scaling relations imply that most starburst galaxies are 
not fully calorimetric for /ssn V = lO""" erg. Our present anal- 
ysis does not tightly constrain the CR calorimetry scenario 
for ULlRGs. Energy-independent CR transport should also 
be examined as a potentially important energy-loss mecha- 
nism for CR nuclei in st arburst galaxies (see in this context 
iPersic & Rephaelil 120 12h . Our observations suggest that a 
substantial fraction of CR energy escapes into intergalactic 
space for most galaxies since CR nuclei are expected to dom- 
inate the CR energy budget, but that some starburst systems 
such as M82, NGC 253, NGC 1068, and NGC 4945 have sub- 
stantially higher calorimetic efficiencies relative to the Milky 
Way. A detailed discussion of particular starburst galaxies as 
potential CR calorimeters is presented by Lacki et al. (201 1). 

5.2. Contribution to the Isotropic Diffuse Gamma-ray 
Background 

Celestial gamma-ray radiation not resolved into indi- 
vidual sources is often treated as the sum two compo- 
nents: a spatially structured component originating mainly 
from CR interactions in the Milky Way, and a nearly 
uniform all-sky component often called the 'isotropic dif- 
fuse gamma-ray background' (IGRB). Shortly after the dis- 
cove ry of an isotropic component using the OSO-3 satel- 
lite (IClark et al.l |1968|) and earl y spectral characterization 
with SAS-2 dpichtel et al.l Il975h . several authors proposed 
that multitudinous extragalactic sources too faint to be in- 
dividually resolved must contribute to the observed flux 
(IStrong et al.lll976t iLichti et al.lll978h . The spectrum of the 
isotropic component has been subsequently mea sured with 
EGRET jSreekumaret al.lll998l:IStrong e t al."2004t). and most 
recently with the Fermi LAT (lAbdo et a l. 2010h). Note that 
the intensity attributed to the isotropic diffuse component is 
instrument-dependent in the sense that more sensitive instru- 
ments are capable of extracting fainter individual sources. In 
this work, the term IGRB will refer specifically to the most 
recent measurement reported by the fen«/-LAT Collabora- 
tion, consistent with a featureless power law of spectral index 
of 2.41 ±0.05 between 0.2-100 GeV with integral intensity 



above 0.1 GeV of 1.03 ±0.17 x 10"^ ph cm'^ s"' sr"'. 

The origin of the IGRB flux is not yet fully understood, 
in part because the contribution from the most prominent ex- 
tragalactic gamma-ray source class, namely blazars, is well- 
constrained by population synthesis techniques. More than 
one thousand sources at high Galactic latitudes, predomi- 
nantly blazars, ha^ve no w been discovered at GeV energies 
jAckermann et al.ll201 lb . The total intensity of these resolved 
sources averaged over the full sky is 0.44 x 10"^ ph cm"^ s"' 
sr"' . Based on the empirically-determined flux distribution of 
high-latitude sources, blazars with fluxes below the LAT de- 
tection threshold ar e expected to contri bute less than 30% of 
the IGRB intensity ( lAbdo et al.ll2010ih . 

Star-forming galaxies far outnumber AGN in number den- 
sity, but are more challenging to detect in high-energy gamma 
rays due to their comparatively modest luminosities and lack 
of beamed emission. Even with the improved sensitivities 
of contemporary gamma-ray telescopes, a limited number of 
exclusively low-redshift star-forming galaxies are expected 
to be individually detected. Estimates for the collective in- 
tensity of unresolved galaxies consequently rely upon scal- 
ing relations between gamma-ray luminosity and galaxy pro- 
porties (such as SFR and gas content), and studies of the 
ev olving cosmological popu lati on of galaxies. Calcula tions 
by lPavlidou & Fields! (I2002t) andlThompson et alJ (120071) dur- 



ing the EGR ET era, and by 



( 12011b . and IStecker & VentersI (1201 ll) incorporating addi 



Fields et al 



tional contraints from the first year of the Fermi mission 
demonstrated that star-forming galaxies could make a sub- 
stantial contribution to the IGRB, comparable to that of the 
bl azars. 

[Fields etaP (1201 Oh pointed out that uncertainties in the esti- 
mated contribution can arise from a degeneracy between den- 
sity and luminosity evolution of star-forming galaxies over 
cosmic time. In other words, a change in either the density of 
galaxies having a given SFR (density evolution), or a change 
in the SFR of individual galaxies (luminosity evolution) could 
account for variations in total SFR density with redshift. If the 
relationship between gamma-ray luminosity and SFR is non- 
li near, such distinctions are important. 

IStecker & VentersI (1201 ll) proposed an approach combin- 
ing global gamma-ray luminosity scaling relations with mul- 
tiwavelength luminosity functions which describe the abun- 
dances of individual galaxies within different luminosity 
classes over redshift. As an example of this approach, we 
now estimate the contribution of non-AGN dominated star- 
forming galaxies to the IGRB using the relationship be- 
tween gamma-ray luminosity (0.1-100 GeV) and total IR 
(8-1000 i-Lm) luminosity identified in section 14.3! together 
with an IR luminosity function provided by Spitzer observa- 
tions of the VIMOS V LT Deep Survey and GOODS fields 
( lRodighieroetal.1120101) . 

This procedure requires an IR luminosity function for ex- 
clusively non-AGN dominated star-forming galaxies. In ob- 
taining luminosity functions from Spitzer data, Rodighiero 
et al. utilized a catalog of spectral templates to simulta- 
neously determine the redshift and spectral classification of 
each galaxy in their dataset. Objects classified as type- 
1 AGN were subsequently removed from the sample. Al- 
though some AGN likely remain, the fractional contamina- 
tion of IR-bight obscurred AG N is expect ed to be small com- 
pared tojion- AGN galaxies (iBallantyne & Papovich 20071 
lPetridl20Tol) . such that the final IR luminosity functions ac- 
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curately represent the cosmic star formation history. The 
IR luminosity function for non-AGN galaxies published by 
Rodighiero et al. (given in the rest-frame of the galax- 
ies) is consistent wi t h those found from o t her IR surveys 
jSanders et al.' '2003; Chapm an et"an l2005t lie Floc'h et al.l 
l2005- , Huynh et al. 2007: Ma gneUi et alJl20m 



The collective intensity (ph cm ^ s ' sr" GeV"') of un- 
resolved star-forming galaxies at observed photon energy £0 
can be computed using the line-of-sight integral, 



(25) 

where $(L^,z) = cfN/{dVdL^) expresses the number 
density of galaxies per unit luminosity interval, and 
dN/dE (^L^,Eo{l +z)) is the differential photon flux of an in- 
dividual galaxy with integral gamma-ray luminosity at red- 
shift z- The factor d^V /{dzdil) represents the comoving vol- 
ume element per unit redshift and unit solid angle. We trans- 
form the integral over gamma-ray luminosity into an integral 
over IR luminosity using the scaling relation equation |4] and 
fitted parameters from Table |3] 

Given the limited number of star-forming galaxies detected 
in high energy gamma rays, the spectral transition from qui- 
escent to starburst galaxies has not yet been mapped in de- 
tail. Accordingly, two spectral model choices for star-forming 
galaxies in the gamma-ray energy band are used in our calcu- 
lation. First, we adopt a power law spectral model with photon 
index 2.2, characteristic of the LAT-detected starbursts. The 
second spectral model is ba sed on a model of th e global emis- 
sion from the Milky Way ([Strong et al.l 120101) scaled to the 
appropriate corresponding IR luminosity. These two spectral 
models should be viewed as bracketing the expected contri- 
bution since multiple galaxy types with different gamma-ray 
spectral characteristics contribute to the IGRB, e.g. dwarfs, 
quiescent spirals, and starbursts. 

Only the redshift range < z < 2.5 is considered in this 
work because IR luminosity functions are not yet well con- 
strained beyond z ^ 2.5. We explicitly assume that the re- 
lationship between gamma-ray and IR luminosity found for 
galaxies at z < 0.05 holds for galaxies at higher redshifts. Al- 
though we are unable to validate this assertion directly, the 
closely-related R C-IR correlation do es not show signs of evo- 
lution up t o z ^ 2 ([Ivison et alJ2010h . and possibly to redshifts 
z > 4 (Sa rgent et al.ll2010h . Those observations are indicative 
of a consistent relationship between star formation and CRs 
in galaxies over the past 10 Gyr. 

We account for the attenuation of gamma rays by interac- 
tions with the extragalactic background light using the model 
of Franceschini et al. (2008). The contribution of resulting 
cascade emission to the IGRB i s expected to be fa int com- 
pared to the primary component (iMakiya et al.ll20lTl) . 

Estimates for the collective intensities of unresolved star- 
forming galaxies, including both quiescent galaxies and star- 
bursts, in the energy range above 0. 1 GeV are shown in Fig- 
ure|2]relative to the first-year Fermi LAT IGRB measurement. 
Shaded regions denote the range of statistical and systematic 
uncertainties. We use the scaling relation between gamma- 
ray luminosity and total IR luminosity obtained using the EM 
algorithm with the full sample of 69 galaxies, allowing for 1 
standard deviation of variation in the fitted parameters, and in- 
clude log-normal intrinsic scatter of gamma-ray luminosities. 



An additional 30% uncertainty is allocated for the normal- 
izaion of the IR luminosity function. Using the methodology 
described above with the scaled Milky Way spectral model, 
the estimated integral photon intensity of unresolved star- 
forming galaxies with redshifts < z < 2.5 above 0.1 GeV is 
1.2^0^ X 10"^ ph cm~^ s"' sr"'. If the power law model is as- 
sumed, the estimated integral photon intensity is 0.8;!;y 4 x 10~* 
ph cm"^ s"' sr"'. These estimates cover a range of 4-23% of 
the integral IGRB intensity above 0. 1 GeV measured with the 
LAT, 1.03 ±0.17 X 10-5 ph cm'^ s"' sr"' dAbdo et al.ll2010hl) . 
Differential intensities for the two spectral models are pre- 
sented in Table |6l The range of intensities predicted in this 
work is consistent with previous estimates for the intensity of 
unresolved star-forming galaxies shown in Figure |7]for com- 
parison. 

The relative contributions of star-forming galaxies to the 
IGRB according to their redshift and IR luminosity class are 
shown in Figure [8] In both panels, fractional contributions 
are normalized to the estimated total intensity of star-forming 
galaxies with redshifts < z < 2.5. These panels provide im- 
portant checks to the reliability of our total estimate given 
the uncertainty in the relationship between gamma-ray and 
IR luminosity for all galaxy types (e.g. no ULIRGs are de- 
tected at GeV energies), and our incomplete knowledge of the 
IR luminosity function at intermediate redshifts. The black 
dashed curve indicates the total IR luminosity above which 
iRodighiero et al.l (120101) report completeness in the sample of 
galaxies used to derive their published IR luminosity func- 
tions. The luminosity function in the region of phase space 
with IR luminosity above that threshold is well-constrained 
by Spitzer data. Although the IR luminosity function for 
galaxies below this threshold is not directly determined by 
the detections of individual galaxies, the total density of IR 
emission is contrained by the observed extragalactic back- 
ground light. The IR luminositity functions presented by 
IRodighiero et al.l (120101) integrated over IR luminosity and 
redshift are consistent with the cosmic SFR history estimated 
from a nalyses of IR observa tions bv iPerez-Gonzalez et al.l 
(120051) . iLe Floc'h etaP (120051) . and ICaputi et al.1 (120071) and 
that estimated with far-ultraviolet data (iTresse et al I I2OO7I) . 

Although the contribution from star-forming galaxies is 
only considered in the redshift range from < z < 2.5 in this 
work, Figure|8]makes apparent that the contribution from star- 
forming galaxies at redshifts z > 1.5 is diminishing. The par- 
ticular shape of the cumulative contribution curve in the lower 
panel of Figure |8] is determined mainly by the IR luminosity 
function. 

Estimates for the combined intensities of unresolved 
blazars and star-forming galaxies are plotted in Figure The 
summed contribution falls short of explaining the IGRB, sug- 
gesting that important aspects of these populations are not in- 
cluded in current models and/or that other high-energy source 
classes or diffuse processes con tribute a signifi cant fraction 
of the observed IGRB intensity. iDermeij (120071) provides an 
extensive discussion of possible contributions to the IGRB. 

Besides the beamed emission from blazars, in which the 
relativistic jet direction coincides with our line of sight, the 
gamma-ray emission from the cores of misaligned AGN must 
also constitute part of the IGRB. More than 10 such radio 
galaxies have now b een detected in Fer mi LAT data at red- 
shifts up to z 0.7 dAbdo et alj|2010el) . A relationship be- 
tween the gamma-ray and radio emission of LAT-detected 
misaligned AGN has been used to estimate that gamma-ray- 
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loud radio galaxies a ccount for ^ 25% of the unresolved 
IGRB above 0.1 GeV ( IInouell20Tlb . However, estimates for 
this contribution are not yet tightly constrained due to the 
small sample of LAT-detected objects and uncertainties re- 
garding the scaling relation between radio and gamma-ray lu- 
minosities. 

Large scale structure formation shocks leading to the as- 
sembly of galaxy clusters represent another interesting can- 
didat e population to explain the remaining IGRB inten- 
sity (IColafrancesco & Blasilll998h . The tenuous intergalac- 
tic medium of galaxy clusters is thought to be a reser- 
voir for CR nuclei which accumulate over cosmological 
timescales. Although non-thermal leptonic populations are 
well-established by observation of Mpc-scale radio features, 
no g alaxy cluster has been de tected in the GeV energy range 
(e.g. lAckermann et al.l 120101) . Other extragalactic source 
classes have been considered, such as gamma-ray bursts 
jLe & Dermeill2007b, and extended emis si on from the lobes 
of ra dio galaxies (IStawarz et al.l 120061: iMassaro & Aiellol 
1201 lb . but these contributions are exprected to be less than 
1% in the GeV energy range. Unresolved radio-quiet AGN 
are also not expected to make a significant contribution to the 
IGRB (Teng et al. 20l3- 

An example of a truly diffuse component is the cumulative 
emission resulting from electromagnetic cascades of ultra- 
high energy CR s interacting with cosmic mi crowave back- 
ground photons (!Berezinskii & Smirnovll 19751) and very-high 
energ y photons in teracting with the extragalactic background 
light jCoppi &~Aharonianlll997l) . The gamma rays produced 
by those cascades are expected to have a relatively hard spec- 
trum. Already, the LAT IGRB measurement has been used 
to constrain this component and to predict the cosmogenic 
ultra-high energy neutrino flux originating from charged pi on 
decays of the ultra-high energy CR interactions (lAhlers et alj 
11)10; Berezinsky et al. 201 1; Wang et al. 20fll). 

Galactic sources, such as a population of unresolved 
millisecond pulsars at high Galactic latitudes, could be- 
come confused with isotropic d iffuse emission as argued by 
iFaucher-Giguere & Loebl (120101) . Part of the IGRB may also 
come from our Solar System as a result of CR interactio ns 
with debris of the Oort Cloud ("M oskalenko & Portei]|2009l) . 

Finally, a portion of the IGRB may originate from 'new 
physics' processes involving, for in stance, the annihilation 
or decay of dark matter particles (iBergstrom et al.l 1200 1[ 
lUlUo et al.ll200llTavlor & Silkll2003l) . 

Studies of anisotropics in the IGRB intensity on small 
angular scales provide another approach to identify IGRB 
constituent source populations jSiegal-Gaskin si 120081) . The 
fluctuation angular power contributed by unresolved star- 
forming galaxies is expected to be small compared to other 
source classes because star-forming galaxies have the high- 
est spatial density among confirmed extra galactic gamma-ra y 
emitters, but are individually faint (lAndo & Pavli dou 20091). 
Unresolved star-forming galaxies could in principle explain 
the entire I GRB intensity without e xceeding the measured 
anisotropy (lAckermann et al.l 120 12al) . By contrast, the frac- 
tional contributions of unresolved blazars and millisecond 
pulsars to the IGRB intensity are constrained to be less than 
~ 20% and ^ 2%, respectively, due to larger angular power 
expected for those source classes. 

6. GALAXY DETECTION OUTLOOK FOR THE FERMI LAT 

The scaling relations obtained in Section l43] allow straight- 
forward predictions for the next star-forming galaxies which 



could be detected by the LAT. We use the relationship be- 
tween gamma-ray luminosity and total IR luminosity to select 
the most promising targets over a 10-year Fermi mission. 

We begin by creating an IR flux-limited sample of 
galaxies from the IRAS Revised Bright Galaxies sample 
dSanders et al.ll2003l) by selecting all the galaxies with 60 /^tm 
flux density greater than 10 Jy (248 galaxies). Next, 0.1-100 
GeV gamma-ray fluxes of the galaxies are estimated using 
the scaling relation between gamma-ray luminosity and to- 
tal IR luminosity. Intrinsic dispersion in the scaling relation 
is addressed by creating a distribution of predicted gamma- 
ray fluxes for each galaxy, assuming the log-normal intrin- 
sic scatter fitted with the EM algorithm for the full sample of 
examined galaxies. Finally, we estimate the LAT flux sensi- 
tivity at the location of each galaxy given the Galactic fore- 
ground and isotropic diffuse background, and extrapolating 
the current LAT observation profile. For each realization of 
the population, we count galaxies with predicted fluxes above 
the LAT sensitivity thresholds at their respective positions 
as 'detected.' The galaxies are modeled as spatially unre- 
solved sources having power law spectral forms with photon 
index equal to 2.2. These modeling choices are focused on 
the search for additional starburst galaxies beyond the Local 
Group. 

Detection probabilities over a 10-year Fermi mission are 
plotted in Figure [TO] for the best candidate galaxies accord- 
ing to the procedure outlined above. The predictions pass the 
consistency check that the best candidate galaxies match those 
which have been LAT-detected. Galaxies with the next high- 
est probabilities of LAT-detection in the coming years include 
M33, M83, NGC 3690, NGC 2146, and Arp 220. This list 
su bstantially overlaps with the top candidates recently named 
bv lLackiet a"n(l2011l) . 

Formally, the galaxy NGC 5128 (Centaurus A) would have 
entered our list of top condidates based on its total IR flux, and 
indeed, NGC 5 128 has been detected at both GeV and TeV en- 
ergies (lAharonian et al.ll2009T: lAbdo et al.]|2010fl) . However, 
the signal is thought to originate mainly from the relativistic 
jet associated with the AGN. NGC 5128 is consequently ex- 
cluded from our list of candidate galaxies. 

We find moderately significant excesses above backgrounds 
at the locations of NGC 2146 and M83 (see Section g]). If 
those excesses are indeed assoc i ated w ith the galaxies (see 
in this context iLenain & Waited (1201 ll) for M83), then both 
might be securely detected after ~ 6 years of LAT observa- 
tions. 

Figure [TT] shows the cumulative number of galaxies antici- 
pated to be detected as a function of increasing mission time. 
The shaded confidence regions are obtained directly from the 
distribution of total galaxies above the LAT sensitivity thresh- 
old as predicted from the multiple realizations of the galaxy 
population described above. T he actual numbers of external 
galaxies reported i n the IFGL (lAbdo et al.ll2010cl) and 2FGL 
dNolan et al.ll2012l) catalogs are consistent with expectations. 
If the simple scaling relationship between gamma-ray lumi- 
nosity and total IR luminosity found for our present sample 
is representative of the larger population of low-redshift star- 
forming galaxies, then we can expect about 10 external galax- 
ies to be gamma-ray detected during a 10-year Fermi mis- 
sion. 

7. CONCLUSIONS 

We examined 64 galaxies selected for their abundant dense 
molecular gas using three years of Fermi LAT data in the 
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0.1-100 GeV energy range. Those results are combined with 
previous studies of 5 Local Group galaxies in a multiwave- 
length analysis incorporating both measured fluxes for the 
LAT-detected galaxies and gamma-ray flux upper limits. We 
find further, though not yet conclusive (P-values ^ 0.05 in the 
most conservative case), evidence for a quasilinear scaling re- 
lation between gamma-ray luminosity and star-formation rate 
as estimated by radio continuum luminosity or total IR lumi- 
nosity. Since contemporaneous star-formation rates are sen- 
sitive to short-lived massive stars, and because gamma rays 
are the most direct probe of CR hadrons in the ISM of ex- 
ternal galaxies, the scaling relationship strengthens the con- 
nection between CR energy content and massive stars which 
ultimately explode as supernovae. 

In the paradigm that SNRs channel approximately 10% of 
their mechanical energy into CR nuclei with kinetic energies 
> 1 GeV, the normalization of the observed scaling relation- 
ship between gamma-ray luminosity and SFR implies that 
starburst galaxies (SFR ^ lOM© yr"') have an average calori- 
metric efficiency for CR nuclei of 30-50% if the gamma-ray 
emission i s do minated by neutral pion decay. As discussed 
in Section [STl this constraint depends upon the assumed en- 
ergetics of SNRs and the stellar initial mass function. Mean- 
while, the hard gamma-ray spectra of starburst galaxies (F = 
2.2-2.3) give preference for energy-independent loss mecha- 
nisms for CR nuclei in starbursts, as opposed to the diffusive 
losses which likely shape the observed gamma-ray spectra of 
quiescent Local Group galaxies. Both hadronic interactions 
and advective transport of CR nuclei should be considered 
as potentially important energy loss mechanisms in starburst 
galaxies. The overall normalization of the relation between 
gamma-ray luminosity and SFR supports the understanding 
that the majority of CR energy in most galaxies eventually 
escapes into intergalactic space. 

The relationship between gamma-ray luminosity and total 
IR luminosity can be used as a low-redshift anchor to esti- 
mate the contribution of star-forming galaxies to the IGRB. 
Using the gamma-ray luminosity scaling relations presented 
here in conjunction with IR luminosity functions, we estimate 
that star-forming galaxies with redshifts Q < z < 2.5 have a 
sky-averaged intensity of 0.4-2.4 xlO~^ ph cm"- s"' sr"' in 
the energy range 0.1-100 GeV, thereby contributing 4-23% 
of the IGRB intensity measured by the LAT. The combined 
contributions of unresolved blazars and star-forming galaxies 
appear to fall short of fully explaining the observed IGRB, 



suggesting that other gamma-ray source populations and/or 
truly diffuse processes constitute a substantial fraction of the 
observed IGRB intensity. 

Finally, we predict that several more external galaxies 
might be detected by the LAT during a 10-year Fermi mis- 
sion. In particular, the galaxies M33, M83, NGC 3690, NGC 
2146, and Arp 220 are interesting targets for further study. 
The all-sky coverage provided by the LAT can also help iden- 
tify promising targets for studies with imaging air-Cherenkov 
telescopes, including the proposed Cherenkov Telescope Ar- 
ra>Q 
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APPENDIX 
KENDALL r CORRELATION STATISTIC 

The Kendall r correlation statistic is a non-parametric rank-correlation method which can b e generalized for the anal ysis of 
datasets containing both detections and non-detections. In this work, we follow the procedure of lAkritas & Siebe"rtl(^1996^ to test 
the degree of correlation between luminosities in two wavebands for a collection of galaxies. 

Consider a dataset consisting of n points, Xk{i], indexed by / = !,...,«, with k = 1,2 to denote the two wavebands. To compute 
the correlation coefficient, first define a helper function, 

Mi J) = 5{Xum{Xk[i\ < Xdj])-S(Xdj])I(Xdj] < XkUm (Al) 

where / = 1 if the enclosed conditional statment is true, and 7 = otherwise. When performing an analysis of partially censored 
data, let S{Xit[i]) = 1 for a detection (measured value), and 5{Xic[i]) = for a non-detection (limiting value). For historical reasons, 
the Kendall r statistic is conventionally defined in the context of data containing lower limits rather than upper limits as is 
common in astrophysics research applications. The method can be applied to data containing upper limits by taking the negative 
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value of all quantities, i.e. = -Xi:[i]. Next, define a function to encode the ordering between two points in the dataset labeled 

H(i,j) = JdiJ)J2(iJ)- (A2) 

A graphical representation of this encoding is shown in Figure [T2]for a dataset containing upper limits in the k = 2 waveband. 
Notice that H{i,j) € -1,0, 1. The Kendall r correlation coefficient is the sum of H{i,j) for each pair of points in the dataset, 
normalized by the total number of pairs, 

2 " " 

T = — (A3) 
n(«- 1) ^-^ ^-^ 

i j>i 

POTENTIAL HAZARD OF MULTIWAVELENGTH COMPARISONS IN FLUX-SPACE 

Consider the simple case of a sample of objects observed in two wavebands, denoted by k= 1,2, which follow a power law 
relation in intrinsic luminosity, log L2 = a logLi + /3. An algebraic conversion of this equation from luminosity-space to flux-space 
yields, 

\0gF2 = alogFi+2{a- l)log(i -I- constant, (Bl) 

where d is the distance to the distance to each object in the sample, = AncfiFk. Notice that if the relationship between intrinsic 
luminosities is non-linear, a ^ 1, the distance to each galaxy enters as an additional term in the flux-space relation. This situation 
can lead to unpredictable behavior, since the objects in the sample are, in general, located at different distances to the observer. 

SIMULATIONS OF MULTIWAVELENGTH COMPARISONS IN LUMINOSITY-SPACE 

We discuss results from a Monte Carlo study of multiwavelength comparisons in luminosity-space in the context of the present 
work. One concern of comparing luminosities is that both quantities have a shared dependence on the measured distances to 
the objects. Common selection effects in observational astrophysics, e.g. that more luminous objects can be detected at greater 
distances, therefore have the potential to induce spurious apparent luminosity correlations between wavebands. 

We performed a series of simulations of a population of objects observed in two wavebands to investigate this potential source 
of bias. The two wavebands are differentiated in that the 'selection' band, here called the 'X-band,' is used to select a sample of 
objects from a larger population, and the objects are then observed in a second 'search' band, here called the 'F-band.' The bands 
are asymmetric in that partial censoring is possible in the Y band, whereas all objects are detected in the X band. The galaxies 
examined in this work fit the above description as a flux-limited sample selected according to their IR, CO-line, and HCN-line 
fluxes, and subsequently observed in the gamma-ray band. 

We test various scenarios considering different intrinsic relationships between X- and K-band luminosities, sample selection 
approaches, and detection efficiencies. Throughout the simulations, we assume that distances to the objects and fluxes are 
measured without error This assumption implies completeness down to the detection threshold flux level in the Y band. 

For each scenario, repeat the following procedure many times: 

1 . Create a population of objects with randomly distributed distances from the observer and luminosities in the X band. 

2. Assume a relationship between intrinsic luminosities in the X and Y bands, specifically a power law form with intrinsic 
dispersion V: 

\ogLY = a\ogLx + l3 + V. (CI) 

Intrinsic dispersion in the relationship between luminosities in the X- and F-bands is taken to be normally distributed in 
log-space with standard deviation axi, i.e. V = M{Q, crp). 

3. Select a subset of objects from the full population. The subset will comprise the sample of objects observed in both 
wavebands and considered in subsequent statistical analysis. 

4. Objects are considered detected (non-detected) in the Y band if their F-band flux is above (below) the F-band detection 
threshold flux level. For the objects which are not detected in the Y band, set the flux upper limit at the F-band detection 
threshold flux level. Determine corresponding luminosities (possibly limits) for each object in the sample. 

5. Compute the Kendall t correlation coefficient for the relationship between luminosities in the X and F bands for the subset 
of selected objects. 

For the particular examples which follow, we generate a population of objects representative of low-redshift star-forming 
galaxies. The X-band luminosity distribution is drawn from the total IR luminosity function of Rodighiero et al. (2010) at 
redshift z = 0, with a minimum IR luminosity of IQ^Lq. The volume considered is a cube of 200 Mpc on a side. Objects are 
assumed to be uniformly distributed in Euclidean three-dimensional space. 
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We consider two sample selection approaches: volume-Umited, and flux-limited in the X band. In either case, a sample 
consisting of 100 objects is selected for each trial. The selected sample in each trial amounts to approximately 0.06% of objects 
in the full simulated population. 

The parameters of each simulation scenario include the intrinsic power law scaling index between X- and K-band luminosities 
(a), the standard deviation of the intrinsic dispersion in the luminosity-luminosity relationship (cr, dex), and the detection effi- 
ciency of sample objects in the Y band. For each scenario, we perform 1000 trials. The distributions of Kendall r correlation 
coefficients for the different scenarios are summarized in Table|2]for the case of volume-limited sample selection, and in Table[8] 
for the case of flux-limited sample selection. 

Volume-limited samples are generally considered to be unbiased, provided that a large enough volume can be sensitively probed 
so as to contain a distribution of objects representative of the full population. The results contained in Table|2]demonstrate that in 
scenarios with no intrinsic relationship between luminosities (a = 0), the distribution of correlation coefficients is consistent with 
being centered on zero (r = corresponds to no correlation). The correlation test is not biased towards finding a non-existent 
correlation in this example. By contrast, the distribution of correlation coefficients for scenarios in which an intrinsic linear 
correlation was considered (a = 1) consistently tend toward positive values, thereby correctly identifying a positive correlation. 

These conclusions do not change appreciably in flux-limited sample selection case, as evident in Table [8] For the simple 
examples above, which are modelled on the population of star-forming galaxies examined in this work in terms of luminosity 
distribution, sample size, and sample selection, the Kendall t correlation coefficient proves to be effective in distinguishing 
between populations with and without intrinsically correlated luminosities, even when only a small fraction of objects are actually 
detected in the 'search' band. 

We considered the simple case of power law relation between intrinsic luminosities for this study. However, since the Kendall 
T correlation test is a non-parametric method, the conclusions are expected to hold for more complex scaling relations as well. 
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TABLE 1 

Summary of the star-forming galaxy sample: beyond the Local Group 
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Galax y distances, total IR (8-1000 /j.m ) lum inosities, and HCN line luminosities are provided by Gao & Solomon i'2004b|). RC luminosities at 1.4 GHz come p rimarily fromf^n et alj 
f200l|) . except for M82 and NGC 3627 (Condon et al. 1990). NGC 4945 1 Wright & Otrupcek 1990), andArp299, NGC 5775, NGC 7331, and VII Zw 31 1 Condon et alj 19981) . Galaxies 
appearing in the Swift BAT 58-month survey catalog with AGN-type classification are identified jBaumgartner et al.|207o|) . 
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TABLE 2 

Summary of star-forming galaxy sample: Local Group galaxies 



Galaxy 


D 
(Mpc) 


LlA GHz 

(10-° WHz"') 


Li2 fim 


L25 fim 

(10^2 


WHz"') 


Lioo 


(10' Lq) 


Lhcn 
(lO' Kkms"' pc-) 


i-0. 1-100 GeV 

(10^* erg s-') 


SMC 


0.06 


0.19±0.03 


0.0029 


0.012 


0.29 


0.65 


0.07±0.01 




0.11 ±0.03 


lmc 


0.05 


1.3±0.1 


0.083 


0.23 


2.5 


2.5 


0.7±0.1 




0.47±0.05 


M33 


0.85 


2.8±0.1 


0.28 


0.34 


3.5 


11 


1.2±0.2 




<3.5 


M31 


0.78 


6.3±0.3 


1.2 


0.80 


4.0 


22 


2.4±0.4 




4.6±1.0 


Milky Way 




19±6 


12 


7.2 


16 


46 


14±7 


4±2 


8.2±2.4 



Global RC, IR, and gamma-ray luminosities of the Milky Way have been estimated using a numerical model of CR propagation and interactions in the ISM iStrong et alJ!2Ql(]|). 
Within the Local Gr oup, an estimate for the global HCN line luminosity is only available for the Milky Way i Solomon et al. 1992). Radio data for other Local Group galaxies: SMC 
iLoiseau et aUl987|) . LMC (Hughe ^^l^200 "^). M3 1 and M 33 ^De nnison et d..l975.) . IR data for other Local Gro up galaxies come from lSanders et al'1i20Q3|) . Gamma-ray data for 
other Local Group galaxies: SMC lAbdo et aU2QlQbl) . LMC jAbdo et alJ20102l) . M31 and M33 jAbdo et alJ2"QlM . 
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TABLE 3 

Maximum likelihood analysis results 



Galaxy 


D 
(Mpc) 


^b.l-lOOGeV 

(10"' ph cm"' s"') 


r 


i'O.l-lOOGeV 

(lO™ erg s"') 


NGC 253 


2.5 


12.6 ± 2.0 


2.2 ±0.1 


0.6 ± 0.2 


M82 


3.4 


15.4 ± 1.9 


2.2 ±0.1 


1.5 ±0.3 


IC 342 


3.7 


<2.7 


2.2 


0.3 


NGC 4945 


3.7 


8.5 ± 2.8 


2.1 ± 0.2 


1.2 ±0.4 


M83 


3.7 


<7.0 


2.2 


0.8 


NGC 4826 


4.7 


<2.9 


2.2 


0.6 


NGC 6946 


5.5 


< 1.6 


2.2 


0.4 


NGC 2903 


6.2 


<2.5 


2.2 


0.8 


NGC 5055 


7.3 


<2.1 


2.2 


1.0 


NGC 3628 


7.6 


<3.1 


2.2 


1.6 


NGC 3627 


7.6 


<3.5 


2.2 


1.7 


NGC 4631 


8.1 


< 1.6 


2.2 


0.9 


NGC 4414 


9.3 


<3.1 


2.2 


2.3 


M51 


9.6 


<3.2 


2.2 


2.6 


NGC 891 


10.3 


<4.2 


2.2 


3.8 


NGC 3556 


10.6 


<2.3 


2.2 


2.2 


NGC 3893 


13.9 


<2.9 


2.2 


4.8 


NGC 660 


14.0 


<3.0 


2.2 


5.0 


NGC 5005 


14.0 


<3.4 


2.2 


5.7 


NGC 1055 


14.8 


<2.9 


2.2 


5.5 


NGC 7331 


15.0 


< 1.7 


2.2 


3.2 


NGC 2146 


15.2 


<6.7 


2.2 


13.2 


NGC 3079 


16.2 


<2.2 


2.2 


5.0 


NGC 1068 


16.7 


6.4 ± 2.0 


2.2 ± 0.2 


15.4 ±6.1 


NGC 4030 


17.1 


<3.0 


2.2 


7.6 


NGC 4041 


18.0 


<3.7 


2.2 


10.3 


NGC 1365 


20.8 


<2.5 


2.2 


9.4 


NGC 1022 


21.1 


<2.2 


2.2 


8.5 


NGC 5775 


21.3 


< 1.6 


2.2 


6.4 


NGC 5713 


24.0 


< 1.8 


2.2 


9.0 


NGC 5678 


27.8 


<2.8 


2.2 


18.3 


NGC 520 


31.1 


<2.0 


2.2 


16.7 


NGC 7479 


35.2 


<6.1 


2.2 


65.2 


NGC 1530 


35.4 


<2.8 


2.2 


29.7 


NGC 2276 


35.5 


< 1.4 


2.2 


15.5 


NGC 3147 


39.5 


< 1.8 


2.2 


23.5 


Aip 299 


43.0 


<3.1 


2.2 


49.3 


IC 5179 


46.2 


<0.9 


2.2 


17.3 


NGC 5135 


51.7 


< 1.5 


2.2 


34.2 


NGC 6701 


56.8 


<2.4 


2.2 


65.2 


NGC 7771 


60.4 


<2.0 


2.2 


62.4 


NGC 1614 


63.2 


<2.1 


2.2 


72.2 


NGC 7130 


65.0 


< 1.3 


2.2 


47.5 


NGC 7469 


67.5 


<2.4 


2.2 


94.9 


IRAS 18293-3413 


72.1 


<2.0 


2.2 


90.8 


Aip 220 


74.7 


<4.4 


2.2 


209.6 


Mrk331 


75.3 


< 1.5 


2.2 


71.9 


NGC 828 


75.4 


<2.7 


2.2 


129.9 


IC 1623 


81.7 


< 1.8 


2.2 


100.6 


Arp 193 


92.7 


<2.6 


2.2 


189.2 


NGC 6240 


98.1 


<2.3 


2.2 


186.5 


NGC 1144 


117.3 


< 1.2 


2.2 


141.1 


Mi-k 1027 


123.5 


<5.6 


2.2 


731.8 


NGC 695 


133.5 


<4.2 


2.2 


646.1 


Arp 148 


143.3 


<5.2 


2.2 


923.8 


Mrk 273 


152.2 


<2.9 


2.2 


582.0 


UGC 05101 


160.2 


< 1.3 


2.2 


287.4 


Aip 55 


162.7 


< 1.7 


2.2 


380.1 


Mrk 231 


170.3 


< 1.9 


2.2 


468.4 


IRAS 05189-2524 


170.3 


< 1.6 


2.2 


395.6 


IRAS 17208-0014 


173.1 


<5.6 


2.2 


1434.0 


IRAS I0566-^2448 


173.3 


<2.0 


2.2 


521.0 


VIIZw31 


223.4 


< 1.6 


2.2 


692.6 


IRAS 23365-f3604 


266.1 


< 1.7 


2.2 


1059.0 



109.4 
180.1 



33.: 



38.1 



Each galaxy is analyzed using a power law spectral model, dN/dE E ^ . Flux upper limits for galaxies not significantly detected in LAT data are presented at the 95% confidence 
level, assuming a photon index r=2.2. 
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Fig. 1. — Gamma-ray light curves for four significantly detected galaxies. The full three-year observation period was divided into 12 time intervals of ~ 90 
days each. A maximum likelihood fit was performed for each of the shorter time intervals to test for variabihty. The dashed black line shows the maximum 
likelihood flux level obtained for the full three-year observation period. The gray band represents a 2% systematic uncertainty in source exposure resulting from 
small inaccuracies in the dependence of the instrument response on the source viewing angle, coupled with changes in the observing profile as the orbit of the 
spacecraft precesses. Flux upper limits are shown at the 95% confidence level. 



TABLE 4 

Scaling relationships between global luminosities: non-parametric analysis of correlation significance 







Full Sample 






Excluding AGN 






(TD = 0% 


O-D = 10% 


CTo = 20% 


(TD = 0% 


(TD = 10% (TD = 20% 


LlAGHi '■ io. 1-100 GcV 


0.001 


0.005 


0.01 


0.007 


0.03 


0.06 


i'g-lOOOpm '■ ^O.l-lOOGcV 


0.0003 


0.001 


0.005 


0.0004 


0.004 


0.02 


i-HCN : ^.1-lOOGeV 


0.05 


0.07 


0.1 


0.2 


0.3 


0.4 



Results for the Kendall r significance tests for correlation between gamma-ray luminosity and each of RC luminosity, IR luminosity, and HCN line luminosity, 
expressed as P- values representing the probabilities of erroneously rej ecting the null hypothesis that no correl ation exists between wavebands. A description of 
the Kendall r statistic is provided in Appendix |A] and details of the implementation are provided in Section 1431 Scatter in the distance measurements to the 
galaxies are assumed to be noraially distributed with standard deviation ((td) expressed as a percentage of the actual galaxy distances. 
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Fig. 2. — Spectral energy distributions for four starburst galaxies significantly detected in high energy gamma rays. 95% confidence upper limits are indicated 
in energy bins with detection significance TS < i. The best-fit power law spectral model is shown for each galaxy in the energy range used for the maximum 
likelihood analysis (0.1-100 GeV; solid) and in extrapolation to neighboring wavebands (dashed). Flux measurements produced by imaging air-Cherenkov 
telescopes are includ ed for M82 (VER ITAS. Acciari et al. 2009) and NGC 253 (H.E.S.S. .Iphm et al. 201 1]), and flux upper limits above 210 GeV are shown for 
NGC 1068 (H.E.S.S.. IAharonian et al.>2Q05.V The uncertainties depicted for LAT flux measurements represent combined statistical and systematic uncertainties. 



TABLE 5 

Scaling relationships between global luminosities: power-law fits with full sample 





Expectation-Maximization Method 


Buckley-James Method 




a l3 V Variance 


a 




vVariance 


Full sample of galaxies 


iL4GHz : io. 1-100 GeV 
(10-' WHz-') : (ergs-') 


1.10 ± 0.05 38.82 ±0.06 0.17 


1.10 ±0.06 


38.81 


0.20 


ig-lOOO iim '■ i«.I-100GcV 
(1O'«L0): (ergs-1) 


1.17 ±0.07 39.28 ±0.08 0.24 


1.18 ± 0.10 


39.31 


0.31 


Excluding galaxies hosting SM>ift-BAT detected AGN 


LlAGHz '■ io.l-lOOGeV 
(10^' WHz-') : (ergs-') 


1.10 ±0.07 38.81 ±0.07 0.19 


1.09 ±0.11 


38.80 


0.24 


is-lOOOjim : io.l-lOOGcV 
(lO'f Lq): (ergs-1) 


1.09 ± 0.10 39.19 ±0.10 0.25 


1.10 ±0.14 


39.22 


0.33 



Fitted parameters for relationships between gamma-ray luminosity and multiwavelength tracers of star-formation. Using the RC case as an example, the scahng 
relations are of the form logLo. i-iooGcV = c«logii.4GHz + /9> with luminosities expressed in the units provided in the leftmost column. The square root of the 
variance provides an estimate of the intrinsic dispersion of gamma-ray luminosity residuals in log-space about the best-fit regression line. For the EM algorithm, 
the intrinsic residuals about the best-fit line are assu med t o be normally distributed in log-space. The Buckley- James algorithm uses the Kaplan-Meier method 
to estimate the distribution of residuals. See Section |43] for a description of the fitting methods. Both the complete sample of 69 galaxies (containing 8 LAT 
sources) and a subsample of 60 galaxies excluding the AGN detected by the Swift BAT are analyzed (containing 6 LAT sources). 
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Fig. 3. — Top Panel: Gamma-ray luminosity (0.1-100 GeV) versus RC luminosity at 1.4 GHz. Galaxies significantly detected by the LAT are indicated with 
filled symbols whereas galaxies with gamma-ray flux upper limits (95% confidence level) are marked with open symbols. Galaxies hosting Swift-BPiT AGN are 
shown with square markers. RC luminosity uncertainties for the non-detected galaxies are om itted for cl arity, but are typically less than 5% at a fixed distance. 
The upper abscissa indicates SFR estimated from the RC luminosity according to equation |2]fYunj;LaL [200 Ij). The best-fit power law relation obtained using 
the EM algorithm is shown by the red solid line along with the fit uncertainty (darker shaded region), and intrinsic dispersion around the fitted relation (lighter 
shaded region). The dashed red line represents the expected gamma-ray luminosity in the calorimetric limit assuming an average CR luminosity per supernova 
of £sN V = 'O'" Ei'g Section IsTt . Bottom Panel: Ratio of gamma-ray luminosity (0.1-100 GeV) to RC luminosity at 1.4 GHz. 
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Fig. 4. — As Figure[3] but showing gamma-ray luminosity (0.1-100 GeV) versus total IR luminosity (8-1000 ^im). IR luminosity uncertainties for the non- 
detected ealaxies are omitted for clarity, but are typically ~ 0.06 dex. The upper abscissa indicates SFR estimated from the IR luminosity according to equation 
ni >Kennicutl.,1998hi) . 
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Fig. 5. — Graphical representation of the method adopted to estimate significances of multiwavelenth correlations using the Kendall r statistic. Total IR lumi- 
nosity (8-1000 fim) and gamma-ray luminosity (0.1-100 GeV) are compared using the full sample of 69 galaxies in this example. Null hypothesis distributions 
of correlation coefficients assuming independence between wavebands are shown for 3 levels of uncertainty in distance measurements ( 1 standard deviation of 
the scatter). The null hypothesis distributions are generated from lO*" permutations of gamma-ray luminosities among the galaxies, requiring that the resultant 
gamma-ray fluxes could have been measureable given the flux sensitivity threshold of the LAT. The correlation efficient of the actual data is represented as a 
probability density to account for uncertainty in measured gamma-ray fluxes. The gray bands indicate the inner 68% and 95% of the probability density for 
the actual data. The congelation significance is estimated by c omputing the fraction of null hypothesis realizations with coiTelation coefficient larger than that 
obtained for the actual data (i.e. the P-value). See Section l43l and Appendix |A]f or details. 
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Fig. 6. — Gamma-ray luminosity spectra of star-forming galaxies detected by the LAT and imaging air-Cherenkov telescopes. 
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Fig. 7. — Estimated contributi on of unresolved s tar-forming galaxies (both quiescent and starburst) to the isotropic diffuse gamma-ray emission measured 
by the Fermi LAT (black points. lAbdo et al]|201Qhl) . The shaded regions indicate combined statistical and systematic uncertainties in the contributions of the 
respective populations. Two different spectral models are used to estimate the GeV gamma-ray emission from star-forming galaxies: a power law with photon 
index 2.2, and a spectral shape based on a numerical model of the global gamma-ray emission of the Milky Way (Strong et al. 201(|). These two spectral models 
should be viewed as bracketing the expected contribution since multiple star-forming galaxy types contribute, e.g. dwarfs, quiescent spirals, and starbursts. 
We consider only the contribution of st ar-forming galaxies in th e redshift range < z < 2.5. The gamma-ray opacity of the Universe is treated using the 
extragalactic background light model of Franceschini e t"al] j20O8l) . Several previous estimates for the intensity of unresolved star-fonning galaxies ai'e shown 
for comparison. Thompson et al. ( 2007) treated starburst galaxies as calorimeters of CR nuclei. The normalization of the plotted curve depends on the assumed 
acceleation efficiency of SNRs (0.03 in this case). The estimates of [Fields et alj (201(1) and by IMakiva et al.l 1201 If) incorporate results from the first year of 
LAT observations. [Fields et al. (2010) considered the extreme cases of either pure luminosity evolution and pure density evolution of star-fonning galaxies. Two 
recent predictions from .Stecker & Ventersi 620 1 11) are plotted: one assuming a scaling relation between IR- luminosity and gamma-ray luminosity, and one using 
a redshift-evolving Schecter model to relate galaxy gas mass to stellar mass. 
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TABLE 6 

Collective gamma-ray intensity of unresolved star-forming galaxies relative to the IGRB 



Star-forming Galaxies Isotropic Diffuse 

Re-scaled Milky Way Model Power Law Model, T = 2.2 
E E-dN/dE E-dN/dE E'^dN/dE 

(GeV) (10-'* GeV cm^^ s"' sr"') (lO^* GeV cm^^ s"' sr"') (Ur** GeV cm^^ s"' sr"') 



0.10 


4.0-15 


4.9-18 


150 


0.16 


4.9-18 


4.5-16 


120 


0.25 


5.4-20 


4.1-15 


100 


0.40 


5.4-20 


3.7-14 


82 


0.63 


4.9-18 


3.4-12 


68 


1.00 


4.2-15 


3.1-11 


57 


1.6 


3.3-12 


2.8-10 


47 


2.5 


2.5-9.0 


2.6-9.4 


39 


4.0 


1.9-6.7 


2.4-8.6 


32 


6.3 


1.4^.9 


2.1-7.8 


27 


10.0 


1.0-3.6 


1.9-7.1 


22 


16 


0.75-2.7 


1.8-6.5 


18 


25 


0.55-2.0 


1.6-5.8 


15 


40 


0.39-1.4 


1.4-5.1 


12 


63 


0.26-0.94 


1.2^.3 


10 


100 


0.16-0.56 


0.89-3.2 


8.6 


160 


0.082-0.28 


0.59-2.0 




250 


0.036-0.12 


0.34-1.1 




400 


0.022-0.078 


0.19-0.63 





Estimated intensities of unresolved star-fonning galaxies using two different spectral model assumptions are compared to the IGRB spectmm measured by with 
the LAT. The spectrum of the IGRB is consistent with a power law characterized by a photon index F = 2.41 it 0.05 with an integral photon intensity above 0.1 
GeVof {1.03±0.17)x 10"^ ph cm - s ' sr ' jAbdo et aU2010l]) . The entries in the rightmost column follow from this parametrization for the IGRB spectrum. 
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Fig. 8. — Relative contribution of star-forming galaxies to the isotropic diffuse gamma-ray background according to their redshift and total IR luminosity 
(8-1000 fim) normalized to the total contribution in the redshift range 0<z<2.5. Top Panel: Solid contours indicate regions of phase space which contribute an 
increasing fraction of the total energy intensity (GeV cm"^ s"' sr"') from all star-forming galaxies with redshifts < z < 2.5 and IO^'Lq < Lg-iooo fim < 10"Lq. 
Contour levels are placed at 10% intervals. The largest contribution comes from low-redshift Milky Way analogues (Lg-iooo ^^m ~ IO'^Lq) and starburst galaxies 
comparable to M82, NGC 253, and NGC 494 5. The black dashed cu rve indicates the IR luminosity above which the survey used to generate the adopted IR 
luminosity function is believed to be complete jRodighiero et al.120101) . Bottom Panel: Cumulative contribution versus redshift. As above, only the redshift range 
<z< 2.5 is considered. 
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Fig. 9. — As Figure|7] but showing the summed contributions of blazars and star-forming galaxies (this work) to the isotropic diffuse gamma-ray background. 
Two different assumed spectral models for the star-forming galaxies are shown. The estimated contribution of blazars is deiived fr om the distribution of observed 
fluxes for high Galactic latitude sources observed by the LAT, which are believed to be dominated by FSRQs and BL Lac objects lAbdo et aU2010ih . 
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Fig. 10. — LAT detection probabilities for individual galaxies using the scaling relation found in Section [431 Galaxies already detected by the LAT are shown 
by blue curves, while other galaxies are drawn with red curves. 
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Fig. 1 1 . — Total number of st ar-forming galaxies anticipated to be detected during a 10-year Fermi mission. The actual numbers of external galaxies reported 
in the IFGL <Abdo et aU2010cD and 2FGL jNolan et aU2012l) catalogs are marked with star symbols. 
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Fig. 1 2. — Graphical representation of the generalized Kendall r rank correlation test. The r con'elation statistic is proportional to the sum of //-values obtained 
for each pair of points in the dataset. See Appendix lAlf or a complete desciption of the method. 



TABLE 7 

Kendall t Correlation Coefficient Distributions: Volume-limited Sample 



Detection Efficiency in y-band = 0. 1 







'^intrinsic — 


f^iiitrinsic — 1 




=0 


o±o 


0.12 ±0.023 






( ■■■ ) 


(5.1) 




0.3 


-0.00021 ±0.012 


0.11 ± 0.022 






(-0.017) 


(5.3) 


o-D 


=1 


-0.002 ± 0.028 


0.092 ± 0.024 






(-0.07) 


(3.9) 


Detection Efficiency in F-band = 1 






Q^intiinsic — 


f^iiitrinsic — 1 




=0 


o±o 


1 ±0 






( ■■■ ) 


( ■■■ ) 




0.3 


-0.0018 ±0.067 


0.71 ± 0.03 






(-0.026) 


(24) 




=1 


-0.0023 ± 0.068 


0.36 ± 0.058 






(-0.035) 


(6.2) 



Mean Kendall r correlation coefficient and standard deviation of the distribution are supplied for each set of simulation parameters: detection efficiency in the 
Y band, true power law scaling index of intrinsic luminosities (a), and standard deviation of intrinsic dispersion in the luminosity-luminosity relationship (cr-p, 
dex). 1000 realizations were analyzed for each set of simulation parameters. Below the mean and standard deviation of the distribution of correlation coefficients, 
and indicated in parentheses, is the ratio of the mean to the standard deviation. 
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TABLE 8 

Kendall t Correlation Coefficient Distributions; Flux-limited Sample 



Detection Efficiency in F-band = 0. 1 



*-^intiinsic ~ *^intiinsic ~ 1 



o-D=0 ± 0.068 ± 0.017 

( ■■■ ) (4.1) 

crx,=0.3 0.00074 ± 0.0056 0.072 ±0.018 

(0.13) (4.1) 

(Tu=l 0.0027 ± 0.017 0.065 ± 0.021 

(0.16) (3.1) 



Detection Efficiency in y-band = 1 



Q^intiinsic — Ctintiinsic — 1 



o--D=0 ± 1 ± 

(...) ( ... ) 

o-D=0.3 0.00004 ± 0.067 0.73 ± 0.03 

(0.00055) (24) 

cr-p=l 0.0018 ± 0.07 0.38 ±0.057 
(0.025) (6.6) 



As Table|7] but consideiing flux-limited samples of objects in X-band flux. 



